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ABSTRACT 


The present work 1S an attempt to put together the most 
relevant aspects of the engineering problems involving dis- 
tributed parameter systems (D.P.S.'s). Simulation and optimal 
control are explained in detail in Chapters II and III. 

The original contribution of this thesis is given in 
Chapters V and VI, where modal control theory and a gradient 
subroutine that searches for the optimal reference coeffi- 
clients are used. AS a result, it was possible to obtain an 
Output distribution better than the one achievable by the 
known methods. This technique works in situations of strongly 
nonlinear control and compensates the effect of having the 
analyzer and synthesizer approximated by low order matrices. 
It also makes it possible to give higher weight to some zones 
of the output distribution in order to have a better local 
fit. The necessary background for understanding Chapters V 


and VI iS given in Chapter IV. 
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I. GENERALITIES ABOUT THE CONTROL OF DIES@RIBUTED PARAMETER 

Svs ters 
Pee HS TORY 

Distributed parameter systems have existed in Nature 
Since the beginning, and many of these systems were natu- 
rally stable. 

For example, during the period of pre-living beings dra- 
matic changes in temperature originated no less dramatic 
variations of the lithosphere. These earth's convolutions, 
together with many subsidiary effects, tried then to make 
changes in different places until a practical equilibrium was 
Obtained. What was this, other than a distributed parameter 
system? 

With the advent of life in this planet other D.P.S.'s ap- 
peared, such as the nitrogen cycle, in which different forms 
of life receive nitrogen and give it up to other forms which, 
in turn, originate new nitrogen. 

In the beginning of this century much attention was dedi- 
cated to the mathematical description of D.P.S.'s such as 
transmission lines. With the entrance on the scene of modern 
control theory and integrated circuits concepts, some attempts 
were made toward a generalization of those ideas, but practi- 
cally until the 1950's not much had been realized. By that 
time a certain number of systems with time delays could be 
analyzed using what are called conventional control techniques, 


but the scope of these techniques was quite limited. 


ILS: 








In the early 1960's some technical papers were published 
in an attempt to apply the concepts of optimal control to 
D.P.S2's and, from then on, a highly geometric rate of tech- 
nical paper production took place. 

A new kind of approach to the control of D.P.S.'s was 
carried on by Murray-Lasso [Ref. 45] and Gould [Ref. 26] in 
1965, using the generalization of the concepts of modal con- 
trol introduced by Rosenbrock [Ref. 53]. This approach is 
treated in depth in Chapters IV, V and VI. References 62 and 
18 also contain very useful information in this field. 

In 1968 the Int. Journal of Control published an exten- 
Sive bibliography by Wang [Ref. 72], one of the scientists in 
the U.S.A. who has contributed most significantly to the de- 
velopment of the D.P.S.'s theory. In November 1969, the 
Aerospace Research Laboratories published the "Survey of op- 
timal Control of Distributed Parameter Systems," by A. C. 
Robinson [Ref. 52], which contains 261 references and an anal- 
ysis of various ways of implementing the optimal control of 
such systems. Special reference ought to be made to the 
contributions of Lions [Ref. 40] and Butkovskiy [Ref. 5], 
Whach respectively in France and U.S.S.R. established founda- 
tions of a scientific treatment of the optimal control of 
Deeso.'s. 

In 1969 and 1970 numerous technical papers on D.P.S.'s 
were published. Many of them may be found in the "Proceedings 
Seehe Joint AWeomatic Contrel Conference,” “Int. Joummal oe 


Control” and “Simulation.” 
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Almost all the above references, although necessary, are 
very theoretical and not much has been done by the control 
engineer in order to implement those methods or to describe 
them in an easy language. It is a purpose of this thesis to 
moll a little of this existing gap. 

History is not only the study of the past events but also 
the careful analysis of them in order to extrapolate for the 
Pieaure. it seems reasonable, on the bastg@ ok past Hhistoryvere 
make some predictions as to future studies of D.P.S.'s. 

1) From a technological point of view it is a known fact 
that the impact of control science in the improvements 
accomplished in the last two decades constituted a "Sine 
gua non" factor for these trends. At the actual rates 
of development of new D.P.S. s (as Mentioned before the 
integrated circuits are a characteristic example) and 
of the capability to describe and control them, it seems 
that within a few years the control science will be ap- 
plicable to almost every different aspect of the existing 
technology. 

ii) The actual and very important ecological problems, namely 
water and air pollution, if mathematically described will 
assume the form of partial differential equations. The 
consequences of being able to deal with such problems in 
a scientific way can be easily implied. 

iii) It seems realizable that partial differential equations 
may be derived that approximate the behavior of some 


segments of society, at least during certain periods of 
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time. When such equations are found it should be pos- 
Sible to use computer optimization (parameter identifi- 
cation) to determine values or expressions for the 
coefficients. These coefficients are but the histori- 
cal constants the historians have been looking for. 

iv) Man is a complex being, and his free nature leads to 
impulsive behavior at random intervals. The effect of 
such behavior on society is often observable. It is 
reasonable to think that the phenomena of man's behavior 
may in some way follow a Gaussian distribution, as do 
many natural phenomena. If equations can be obtained to 
represent society as a D.P.S., then it may be possible 
to treat the impulSive aspects of man's behavior as 
GausSian noise and perhaps the consequences of such be- 
havior can be minimized by building into the social sys- 
tem something similar to a Kalman filter. 

v) AS progress is made in developing equations for various 
Sub-systems within the social systems, it may well be 
that some sub-systems may be found uncontrollable, others 
GOontrollable. Then the theory of D.P.S."s may welblue 
come an important tool for redesign of social systems. 

Finally, in order to give a geometrical interpretation 
fameeme behavior of a D.P.S., Figs. 1.1 and 1.2 are inciuded, 
which show how each single input contributes to the value of 


the output at every point. 
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B. HOW TO CONTROL DISTRIBUTED PARAMETER SYSTEMS 

The ordinary differential equations used by the control 
enginéer to describe lumped parameter systems, with time as 
independent variable, did not prove sufficiently accurate to 
deal with complex systems where either the input, or the out- 
put, or both, are also function of other parameters, generally 
of spatial nature. Therefore, a new theory was developed in 
order to extend some of the ideas so well understood when 
dealing with ordinary differential scuseiene to the much more 
difficult problems described by partial differential equations 
and integral equations. 

There are four basic approaches to the control problem of 
D.P.S.'s; however, none of them is a general method applicable 
to all types of these systems: 

1. Use of Time-Delay Techniques 

The use of time-delay techniques [Ref. 64] has been 
common when dealing with problems such as those that arise in 
long pneumatic control lines, heat transfer in nuclear plants, 
production delay in assembly lines and in general to every 
kind of control systems where time-delays are involved. The 
conventional analysis and design techniques applied to these 
systems requires too much work and for this reason they are 
now studied almost entirely by computer Simulation. 


2. Reduction of the Partial Differential Equations to 
Ordinary Differential or Difference Equations 


Once the ordinary differential or difference equations 


are obtained lumped parameter techniques must be used. From a 
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control point of view the most representative: methods used to 
obtain the above equations are as follows: 
a. Space Quantization (DSCP) 

The space derivatives are replaced by finite dif- 
ferences but the time derivatives are maintained. One of the 
biggest advantages of this method is that it permits an analog 
modeling of the system. The accuracy can be improved by in- 
creasing the number of sections, which must be adequately 
isolated. When the number of sections becomes too large for 
the required accuracy, another method must be chosen. 

b. Time and Space Quantization (DSDT) 

Use 1S made of the well-known techniques for the 
numerical solution of partial differential equations. Al- 
though quite valuable from an analytic point of view and for 
Simulation, it seems a little lengthy if one is trying to ap- 
Pey it to control design. 

c. Laplace Transform 


The Laplace transform method uses the fact that 


pee, t), = oF (xX,S) OE NEG ee ee a _ 
A Sa = —xx «Cand ae cma = Sloe so) f{x,0 ), when 
eee ft (x,t) is transformable, L{f(x,t)] = F(x,s) and oe leet) 
exists. 


The problems that can be treated by Laplace trans- 
formations are quite numerous, the principal objection being 
the intricate expressions to which one arrives when dealing 
with systems somehow more complex (i.e. multidimensional 


systems). 
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SB. ModaiCcont rel 

The modal control approach, Refs. 26 and 45, consists 
in replacing the partial differential operator with an in- 
finite order matrix which is approximated by another matrix 
of order N<o. Through adequate matrix manipulations it is 
possible to uncouple the interconnected inputs and outputs, 
such that separate loops are established, each of which can be 
treated by lumped parameter control techniques. 

Its most noteworthy advantage is the great insight that 
it can give to a problem. The range of operation is, however, 
limited because it can be applied only to completely continuous, 
bounded and linear or linearized operators.~ Also, sometimes, 
the truncation is problematic, either because some of the 
higher order eigenvalues may be of fundamental importance or 
because the convergence of these eigenvalues may be too slow. 
In this case the system would require a high order matrix 
representation, which is not very acceptable in practice. 

The last drawback seems to be the difficulty that often arises 
in the computation of the eigenvalues and eigenfunctions. A 
detailed analysis of this theory 1S given in Chapter IV. 

fa Optimal Control 

Optimal control is, by far, the technique about which 
the largest amount of work has been published. This comes from 
the fact that the control of systems described by partial dif- 
ferential equations varies remarkably with changes in the 


Wattial and boundary conditions, with the definition of cost 


Se See Appendix B for definitions. 
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function and from system to system. Some more or less gener- 
alized solutions have been proposed, normally according to the 
classification of the system operator within the three fol- 
lowing categories: elliptic, parabolic or hyperbolic partial 
differential operator. Optimal control will be discussed in 

a certain depth in Chapter III; it seems to be the approach 

of greatest future. This is due to the increasing knowledge 
on how to describe the cost function, giving convenient 
weight to each of its components, and also due to the tendency 
for the control systems to become highly computerized because 
@emene low cost of integrated circuitry. 

The four basic approaches just described are not com- 
pletely independent of each other. Effectively, in many 
circumstances, all of them are ultimately reduced to the solu- 
tion of problems involving ordinary differential or difference 
equations. There is, however, an important reason for the 
chosen criterion, which can be better explained giving the 
following example. Consider para. 2; it is obvious that after 
the transformation is accomplished the lumped parameter optimal 
control techniques can be used; the difference from para. 5 is 
in the way the optimal control is derived. In the first case 
the optimal control law is derived from an ordinary differen- 
tial or difference equation, while in the last one it 1s 


derived directly from the given partial differential equation. 
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IT. SIMULATION METHODS FOR DISTRIBUTED PARAMETER SYSTEMS 


A. INTRODUCTION 

The distributed parameter part of a physical system is 
most accurately described by partial differential equations. 
Analysis and design of such systems requires solution of the 
partial differential equations. Solution by analytical meth- 
ods can be quite laborious, and such solutions are not help- 
ful in design problems; indeed they may not be convenient for 
many system analysis problems. 

An alternate approach is the simulation of the distributed 
parameter system in a computer (analog, digital or hybrid) and 
the use of this simulation model in computer studies planned 
for analysis and/or design. A number of different simulation 
techniques have been developed for this purpose. Some of the 
methods have been chosen to fit a specific type of computing 
facility, others have been designed to solve a specific class 
of partial differential equations, and still others are general 
methods that can be applied to a variety of problems. 

The purpose of this chapter is to classify and list most 
of the available simulation methods, giving also some of the 
important facts associated with each of them. No detailed 
comparison or evaluation is attempted since the choice of a 
method is guided largely by available computer facilities, 
the nature of the specific problem to be solved, the accuracy 


needed in the solution, and the constraints on cost and time. 
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However the chapter is concluded with a summarizing table 
which indicates some of the pertinent considerations associ- 
ated with each method. This table could be used as a start 


Sie a critical comparison and evaluatvene 


B. <A GENERALIZED CLASSIFICATION 

HE can be said that, in general, Ehe partial dureeren see 
equation describing a distributed parameter system is part of 
an overall control system. In order to analyze the behavior 
of this complete system it is necessary to know how to solve 
the corresponding p.d.e. This can be done either by analyt- 
ical methods or by simulation. 

Simulation has the great advantage that once the system is 
represented for a specific input and initial conditions any 
change of conditions 1S easily realizable; such is not the 
case for the analytical computation. 

The techniques for dealing with multidimensional systems 
are only extensions of those used for unidimensional ones. 
Due to the fact that the simulation of those systems involves 
an enormous amount of hardware, if implemented by analog meth- 
ods, most multidimensional simulations have been done in the 
digital computer and sometimes in the hybrid computer. They 
are, however, tremendously time consuming and for this reason 
only the two-dimensional systems can be modeled with reason- 
able, if not great accuracy [Ref. 49]. The objective of a 
Simulation will be to have for specific positions the output 


as a function of time or, for each instant of time, the output 


Za 








as a function of the coordinates. In a two-dimensional sys- 
tem the solution is normally obtained keeping two coordinates 
fixed; changing the other one and repeating this for several 
combinations of points and time instants until a complete grid 
is obtained. Actually some oscilloscopes are already capable 
SeeorOviding multidimensional pictures. Wins thesis wat eeeal 
only with one-dimensional systems, for which the following 
methods of solution may be considered: 


DSDT - discrete-space, discrete-time 


DSCT - discrete-space, continuous-time 
CSDT - continuous-space, discrete-time 
TSCT - transformed-space, continuous-time 
TSDT - transformed-space, discrete-time 


Schuchmann [Ref. 59] gives a brief and clear analysis of 
these methods, involving error and stability considerations, 
as well as the relative advantages of each one. He doesn't 
consider the three possibilities resulting from taking the 
time Laplace transform, explaining that this would imply the 
need for getting its inverse; this, although feasible is very 
time consuming. However, when considering the solution at few 
points (sometimes only one) the technique seems to be good, in 
particular if using infinite product expansions [Ref. 22]. 
Such technique gives the solution as a product of terms 
(truncated infinite series) in "s", each of them easily im- 
plementable in an analog computer; it will be discussed in more 
detail in Section F of this chapter. 

1. DSDT - Discrete-Space, Discrete-Time 

Mith the very fast digital computers actuallyve wre. 


tent and also because of their enormous accuracy, popularity 
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and availability, most of the work has the tendency to be 
done by the DSDT method. Also with the actual video-display 
Mmits 1t 1S possible to get CRT pictures@ama change veacmin 
the parameters without need for an analog or hybrid computer. 

The basic concepts of the discretization are very 
Simple and, because of its importance, all of the next sec- 
tion will be devoted to this problem. 

Ze DsCT — Discrete-Space, Continues —11 me 

This is the method generally used when there is avail- 
able only the analog computer and the transformation tech- 
niques described previously in this chapter are difficult to 
implement. If the number of nonlinearities or the required 
accuracy is large, the lack of enough multipliers or other 
components in the computer does not permit Sellinehe. The 
avoid such situation the multiplex method has been tried; it 
consists in switching on and off sequentially, with digital 
peomals, the parallel operation of the circuits corresponding 
to each of the discretized positions and in using the same 
analog hardware to simulate the nonlinearities of each branch 
of the parallel combination. 


In Fig. 2.1 is shown how the heat equation 
2 
Zz Xx 
Ox 


os 
ox 


'@ 
| 
il 


(2315) 


: Some companies work with huge analog computers, as for 
example English Electronics Co. which eecua ioe PIL Sle) = 
amplifier Saturn computer. 
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can be simulated in the analog computer. Cis the heat ca- 
pacity and the heat flow F is given by 


du 


ES ka (2.2) 


If the temperature “u" is measured at integer stations and 


the flux at half-integer stations, Eq. 2.2 may be written in 


uu 
: om Nn 2he : 
discrete form as Pele = ee ee ee and Eq-) 2- lesan 
ia —-F 
8) apo) Ol. Rares GC) a sy 2 iby 
become s(kaz), = ~dxln 7 ie 
from which 
C du, _ Rn) mien ele (2.3) 
eae Ax i 


The above procedure may be easily implemented accord- 
ing to the schematic diagram on the following page. 
a §~CSDT — Continuous-Space, Discrete-Time 
This technigue is advisable for the simulation of the 


flow equation 
— + Vaes = £ isc) 


which, in discrete form, 1S written as 


du_ (x) fe) 


i 
sg = = ee a es) —u, 3 (3) ] Vr 





(2.4) 


From this equation it is obvious Ehat At can be de- 
creased without increasing the number of components, the only 
limitation being the stability requirement (1/vAt<l). The 
analog computer with adequate number of track-hold units or 
the hybrid computer are the most efficient ways of solution of 


this problem. 
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Is coe falnetee 2.1. Analog Computer Simulation of the One 
Dimensional Heat Equation 
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=: .29C? —_Trans formed Space Conan eae 

[SCT 1S a method applicable temsyctemns represented by 
Self-adjoint operators (i.e. heat ancewaye eguations) and, 
through convenient assumptions also to non-self adjoint sys- 
eems>. Lhe concepts involved are intimately related... eamenoce 
of the modal control which will be studied in depth in Chap- 
ters IV and V. The most remarkable example of this method is 
the Bubnov-Galerkin transformation described in Appendix A. 
Another representative example is the one developed and ex- 
plained in very comprehensive terms by Hong and Larson [Ref. 
31] which can be applied to any set of first-order partial 
differential equations in normal form. It gives the solution 
as a continuous function of distance for any value of time. 
The only trade-off is the requirement for the computation of 
the inverse transform. However, because the output is 
U(t,s)=£,-£,s+1/2! (£,s°) - -.-, where £,(n=0,1,2,...) are ex- 
clusively functions of time, the inverse transform is very 
easily obtained. In the example given in the paper the analog 
solution of a single nonlinear first-order p.d.e. required 
only 2 integrators, 3 multipliers, 5 amplifiers and 1 divider. 
Because of the excellent features of this method a more de- 
tailed analysis will be given in Section E. 

5. TSDT - Transformed-Space Discrete-Time 

TSDT is an extension of the above method and it is 
becoming increasingly more popular. One example of such a 
method is implemented in the computer run No. 7 of this 


maesis. 
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C. DIGITAL COMPUTER SOLUTION OF PARTIAL DIFFERENTIAL 
EQUATIONS (DSDT) 


The principal source of information )r@r thus section 
the work, "Numerical Solution of Partial Dafferential Equa-— 
tions," by G. D. Smith [Ref. 63], an excellent book written 
in very comprehensive terms. For sophisticated problems more 
advanced texts are recommended, such as Refs. 1 and 17. Due 
to their relevancy only second-order p.d.e.'s will be 
considered. 


ite Parabolic Equation” 


| | Snes Cee 
Consider the parabolic p.d.e. a = K-—, which 
ox 
. | : oe _U 
normalized using the transformation t=kT/L’, x=; and u=5-) 

0 
where Up is some given value of U (generally maximum or mini- 
mum) at zero time. The resulting equation is 

2 

dU _ du 

Re > a (2.355) 
x 


Written in discrete form the above equation becomes 


| Bo s li. e aa s s e s 
Meee ed A ee (2.6) 
k h2 





: Given the second order quasi-linear p.d.e. 
in Bat 
ax—- + b + c—= + e=0 where a,b,c,e may not be functions of 
3,2 dxd0Y ay* 


Second or higher order derivatives (according to the defini- 
tion of quasi-linearity), if b2-4ac is equal to zero the 
Pel@etion is said to be parabolic, if greater than zero is 
hyperbolic and is elliptic when less than Zero. 


Sil 








oa 
u e e = e ® + e e ona ° r) « 
i,jtl ~ “i,3 eM a aus 5 + Ui 41,5) (Qu) 


Geometrically the interpretation of this formula is Simple, 


as Can be seen from Fig. 2.2. 





Pigure 2.2 Grid Solution or 2 Parabo le sbibske pential 
Equation 


Here, all the values at t=0 and x=0 or 2 are Known. The 
computation is, therefore, straightforward and does not re- 
quire iterations. In the case of mixed boundary conditions 
of the type autbStec, where n is the outward normal to the 
Surface and a,b,c are functions of the space coordinates and 
time, only the values of u at t=0 are known. In this case, 


the solution of the problem becomes only slightly more dif- 


ficult, since it is necessary to evaluate with a simple 
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formula the boundary at x=0 for alll disecretestimes, eae ee 
Start of the computation of the elements of each row (Ref. 
47]. 

This may be done introducing a new difference equation 


based on the boundary condition, as for example 


a Ol Siysieelea 
Ll f = —V f f = 
meas 0,4 2 25x Be 


1S more accurate. 


0-9" which 
The last equation has ae as UnkWewn and there tore 
there 1s need for still another one assuming that the heat 
equation is satisfied at x=0. 
The explicit method, although the simplest one, has the 
inconvenient of only being convergent for 0<k/h*<1/2. For 


this reason some other methods were developed which are un- 


@onditionally convergent for all values of r=k/h*. 


The most popular one is the Crank-Nicolson implicit meth- 


od, which is represented by the equation: 


Ba G41 jd tL iy OME tL “Min, 341 
K a we 


eg dee tala 


nh 


A (2.8) 


Or 


ag) ru, pee ee ee 


ee 1, j+1 se ee eee eee 

(2.9) 
and assuming that the system is discretized in N-l sections, 
for each instant of time, N-1l simultaneous equations must be 
solved. The best way of solving this system of equations is 


using any of the well-known iterative techniques applied to 
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tne digital computer. The best known techniques are Jacobi, 
Gauss-Seidel and successive over-relaxation (S.o.r.). The 
first one is very Slow and for this reason, in general, only 


tne other two are used, rewriting Eq. 2.9 as 
Beorl “i,j ~ prt (uy Lb, jt 174d se ae ee 

+ (W513 - 2u; 5 + Urey A (2 akG®) 
and dropping the index jtl, Eq. 2.10 may now be expressed as 


1 


u, = sr (us _y - 2u, + Us4y) + b. (2. aly) 
where b. = u. . t Solan, > = 2. Seer ue ) (2k) 
a eg 2 dis Ja liye tk, : 


According to the above notations and also denoting the itera- 


tions by superscripts, the Gauss-Seidel scheme is obtained: 


b. 
(gts) a —— (n) 
“al sTirey YY eo Tee ey 
ife converges for all r>0. In matrix form, making p = SHEST 
(n+1) (n) 
uy 0 0° 0 --- i 
(n+1) - (n) 
Uy = 0 0 p oe ai) + C (2.14) 
4 (nt) 0 a ne O 0 ny 
+3 ! ; 
' (n+1) 9 oN-l Ne () 


N= 1 ‘en pees “N-1 


The successive over-relaxation scheme is the fastest one and 
is obtained directly from Eq. 2.13 by adding and subtracting 


u, (n) to the right-hand side. - 


The following expression is obtained: 


b. 


fie) 6(n) a Vor ab) 6m) 
i an { i (iS) 


(n) 
i ue Ceram we sean -u, | (2.15) 


wnere w is the relaxation factor (usually 1<w<2). 
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For maximum rate of convergence 


i — a where wu = = cos 7 (2. 6y) 


an b l+r 
Teac 


2. Hyperbolic Equations 





The most commonly used method for the numerical solu- 
tion of hyperbolic equations is the method of "characteristics" 
and it consists in integrating the given p.d.e. along chosen 
emeections (when dealing with two-dimensional systems; 1 
three dimensional systems are considered the integration 1s 
in chosen planes) for which the partial derivatives are re- 
duced to ordinary derivatives. 

The mathematical explanation is as follows: given 


the quaSi-linear equation 


a ar ar 
+ ea hae et = Ps 
aa oo “32 e 0 ( ) 
2 dU ee! eat an 
with b -4ac>0 and defining ear Pray = Ge aay tt TROY = Ss 
an 
and .—5 = € ,», can now be written: 
fe 
dp = ae ax + SP. ay = rdx + sdy 
and 
dq = 92 gx + 294 ay = sdx + tdy (omic) 
9 oY 


peace from EG. 2.17: artbstctte = 40; which 1s susea slo, 
@eeher with Eq. 2.18 to eliminate r and t, giving the final 


expression 


2 
Giomnrel d d 
scady” - Gd +o} - (AB a Me eM - 0 — 2.as) 


15) 








Choosing SY such that the polynomial inside the left { } 


is zero, it will be necessary to solve only 


dp dy dg dy _ 
asx ase SBE te ees = @) (25.20) 


where so is known and takes two real values. 

A detailed procedure for the computer implementation of 
this technique is given in Ref. 63, pp. 101-102. 

The method of characteristics 1s the most accurate one 
for solving hyperbolic p.d.e.'s and whenever discontinuities 
are involved is the one that must be used. However, if no 
discontinuities are present, a convergent finite,differences 


method is much eaSier to implement and will give generally 


adequate results. One scheme that may be used when dealing 


| Oct San 
with the equation ~~ = — =~ is 
0.2 Zz 
xe ot 
m2 Se 2 2 


with rah, which is convergent for r<l. 

Se Elliptic Equations 

Laplace's and Poisson's equations are the two best 

known equations of the elliptical type. These equations have 
the peculiarity of being always integrable in a closed area 
(if a Ee enoonal equation) or volume (if a three- 
dimensional equation) on which boundaries the function or its 
gradient are known. 


For example Poisson's equation 


2 ee. 
aot 25 = £ (X,Y) (2522) 
Ox oY 
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may be represented in normalized fornmemenna mer, y= 


and ut - One possible discrete-version may be 
ae 


at 


2 
s sar s s e * s e ne s e = 
eee, Cpe wegen Sai au; hv f(xy.) (2.2 3) 


which must be verified for every point inside the boundary. 
In general, if there are N points inside the boundary a set of 
N equations will result, linear or nonlinear, in the same way 
as the original p.d.e. 

When the boundaries are curvilinear and the function val- 
ues or its derivatives cannot be represented accurately by 
the chosen mesh, either a finer mesh must be used or the fol- 


lowing equation: 


5 ey eG ne cl ate. si en : i 22 4 : yO i hf, 
il 1 2 D i 2 1 De: 
(2.24) 


according to Fig. 2.3, where fy 16 the value of =f abw(0,0) 







(O<@<1) 
(0@,<1 ) 


Eeeere 2.3. Grid for an Elliptic p.d.e. with Curvilineat 
Boundaries 


oy 








DB. APPLICATION OF THE CONCEPT OF =CHARAGPARIST (Gs 0 etese ar 
COMEUMIERS DSC ORVeSpn) 


There exists many different ways of solving p.d.e.'s in 
a hybrid computer. Maybe the most representative of these 
procedures is the extension of the concept of the character- 
istics described in the last section. Such an extension was 
Originated by Vichnevetsky and Associates [Refs. 69 and 70]. 
It is explained as follows: given the p.d.-e. 


du ou 
oye £(x,t) >> = g(u,x,t) (2.425) 


with given initial and boundary conditions and making = 


men,t) it is possible to write 


du_ ou, ou, dk _ du, -du 

de (OG Oe ene ieee (eo, 
Therefore 

So = Gu; te) (25527) 


which is an ordinary differential equation, integrable by 
analog computer. 

Bemeethe boundary condition Uu(0,t)) Eq= 2.27 ‘canbe wences 
grated at high speed while the integration of x is done 
slowly (characteristic line). Geometrically, for L the 
length of the system this method is exemplified in Fig. 2.4. 
If the highest integration speed is much greater than the 
lower one the time can be taken as a constant for each fast 
sweep; otherwise a correction can be implemented. For fur- 
ther information a specific case is worked out with great 


Geotail in Ref. 70. 
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Figure 2.4. Geometrical Representation of the Method 
of Characteristics in Hybrid Computer 


HeeeeeNALOG SIMULATION OF UNIDIMENSIONAL D.P.S.'s BY 
LAPLACE TRANSFORMING THE SPACE COORDINATES (TSCT) 


Following the same procedure as used by Hong and Larson 


in Ref. 31, consider the system defined by 


au i‘ £(t) se L gle = & (2.28) 
with initial and boundary conditions 
(BEG) 4) Erni SOs (2.29) 
Wit 0)) = <2 


Normalizing the length and by defining u (t,1)=b the distance 


Laplace transform of u will be 


U (t,s) = MGs (2. 30) 


Oo*— 
0) 
WY) 
* 
c 
Oo, 
* 
lI 
Oo‘ 
0) 
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and, after substituting in Eq. 2.28 it comes out? 


du(t,s) 
dt 


-sx du(t 


ny eee) (ee setS) axj+g (t) Uls) = 0 O51. 


O 
or, because the integral is identical to zero for X greater 


iman 1: 


eo(t,s) 


aE + (g (t) + sf (t)) U(t,s) = £(t) [a-e “b] (Zep 


immorder to obtain the I.C.'s for the solwtion of Eq. 2.32 


make 
©” sx 4 -sx -g(0)x 
We, o) = | e u (o,x) dx = fe e 9 dx 
O O 
J 1 Eee Suey!) 
= Stg(o) {l-e ] (25) 


which expanded in McLaurins's series becomes 


: 2 3 4 
Src) (stg (oO) ) (stq(o) ) Soo) pune 
(2234) 
or 
U(t,o) = hy (g(o))-h, (g(1)) sth, (g(2)) 8” _ ees) 
> Remember that L,[24) = f 78% 84 ax = 2 u(t,s) 
xo t Zs ot elie J 
ou _ —Soeeo lim sees Cbs 
ee ox ! ax OX = Dawe ax OX 
| : P P 
s ee he SX ul +s f e792 Gea 
=>-0o 
O O 
= s | e ?* uy dx - Wc, OF 
O 


su(t,s) ~u(t,0) 
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Equation 2.32 for the derived I.C.'s is known to have a 


Series solution of the form 


ij n 


_o n _ 2 
n=o 
where h(n=0,1,2,...) depends only on time. 


Substituting Eq. 2.36 in Eq. 2232 "and equating seniemeoon= 
ficients of equal powers of s provides one set of ordinary 


differential equations: 


anG 

noe + gh, = f (a-b) 

dh 

qe + ghy — £(h -b) 

dh 

aoe + gh, = £(nh__,-b) (2.37) 


The I.C.‘s for this set are derived setting s=0 in Eq. 2.36 


and solving as indicated for the case of hy (0): 





hg (0) Up 
h, (0) = SS 
a ee 2 2 
2 g S t2g stgo 
(l- <2 + = - ---)-(1-5, - 52+ 
Pan. oy oe 2) oy 
= : - ---) 
(257318) 
The following set is obtained: 
g g” g° 
O O 
hy(0) = l- 57+ 37 > qr tc 
7 oa (‘e C223) 
a O O O 
eho) > > 35° as = Sine 
Hong and Larson prove that AL may be obtained from the 
moment expression ha f x u dx. 
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The set of ntl Eqs. 2.37 has n+2 unknowns and, for this 
reason, another relationship is required to define a unique 
solution. 

The approach followed by the authors of the paper is to 


consider the first moment 


Ps ' i ve 
h, = u dx = xu a i ee OLE 
ne 0 o 
1 du 
= u(t,l) - J Xx a dx 
1 du 
ees ! X ay dx (2246) 


Given the continuous nature of x, Eq. 2.40 may be written, 


according to the mean value theorem as: 


1 . 
_ _ du 2Gaee. 3 
where l 
! xu dx h, 
xX = al = a (2.42) 
fou dx 
0 


The function b is obtained from the above equations and 


it is given by 


= Fea (2.43) 


Equations 2.37, 2.39 and 2.43 are the only requirements 
for an analog solution of the problem. 

The example worked out on the reference paper considers 
only two equations of the set and even so shows good results. 


As stated in Section II-B, the important point on this 
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Simulation is the enormous reduction of hardware it aft forwdse 


This method, if implemented digitally, leads to the TSDT 


Fipproach . 


F. SIMULATION BY INFINITE PRODUCT EXPANSIONS (DSTT) 

This technique is based on the Mittag-Leffler theorem 
Pimech states that a function analytic at the origin and wien 
an infinite number of simple zeros, may be decomposed in one 
ierinite product of simple factors. 

Considering the transcendental function cos z, it is 


equivalent to the product (1-2z2/mT) (1+22/m) (1-22/3m) 


(1+22/31) ----- , OX, if in Taylor's series torm- te theese 
a a" 
quence 1 - ee ir ee The first representation is char- 


acterized by the exact preservation of the zeros. 

Goodson [Ref. 21] works out some examples for different 
mypes Of p.d.e.'s. He considers first the equations in the 
normal form and takes the time Laplace transform in order to 
obtain an ordinary differential equation (with s as parameter) 
in state variable form. Once obtained the general solution 
for the given linear boundary conditions, the infinite product 
expansion iS applied to the transcendental functions in the 
solution. The comparison between the product expansion and 
the eigenvalue or Fourier expansions as indicated in Figs. 
2.5a and 2.5b, allows the statements: 

(i) Both expansions have the same eigenvalues. 


(ii) The product expansion preserves the extremum values 


@i the exact solution (no general proo: Dut 26 aoe 
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Figure 2.5a. Eigenfunction Exeansion 


Ee 
2 
~ 
2 
O 
Y) 
U 
fugue 2.55bD. Iniimite Produces ans1 om 


verified in several cases), while the Fourler expan- 
Sion Minimizes the mean square error. 

The product expansion is very Suitable for analog sim- 
ulation due to the nature of its terms. As an example, 
the solution of the one dimensional heat transfer equa- 
tion with pointwise control at x=0, will yield an 
expression as follows: 


u(l,s) = F(s) If — , where u(1,s) is the 


= 2 
n=1 lts/k, 
Laplace transform of the temperature at x=1 and F(s) 


and KO are obtained from the data. 
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(1v) The poles of the transfer function are preserved in the 
case of the product expansion and, therefore, the tech- 
ree 1S recommendable for the simulation of DPS's 
represented by low-order uni-dimensional (otherwise the 


analytic computations would become too long) p.d.e.'s. 


G. MONTE-CARLO SIMULATION 

Although the study of stochastic models is a little be- 
yond the purpose of this thesis, its actual importance 1S 
Such that some considerations are supposedly worthwhile. 
Two basic papers were published, Refs. 10 and 41, the first 
one describing an analog computer approach and the last one 
a hybrid computer solution. For a good understanding of 
these papers it iS necessary that the reader be familiar with 
the stochastic simulation of systems described by ordinary 
differential equations and also with the theory of Markoff 
processes. An adequate theory of both these fields may be 
found in Refs. 38, 39 and 48. 

The great advantages of this method is that is requires 
only a very limited amount of analog hardware or digital mem- 
ory and that it is only necessary to make the computation at 


the points where the solution is desired. 


PoeeOlLHER SIMULATION TECHNIQUES 

Many other distributed parameter systems have been simu- 
lated in slightly different ways that can, however, be 
classified within the description of the above mentioned pro- 


cedures. Among these special mention is made to Refs. 57 


45 








and 58 which use respectively some known results in the theory 
of Bessel functions and a reformulation of the equations des- 
cribing a parallel-flow heat exchanger process. By consider- 
ing an observer riding a section of the fluid, with the second 
method it was possible to reduce significantly the amount of 
hardware required by the conventional finite differences an- 


meg Simulation. 


Lt. CONCLUSION 

After having described so large a variety of techniques, 
the question arises: Which of them is the best?. 

The answer is not a Simple one. Depending on the avail- 
able computers and also on the capabilities of each one (speed, 
size of the memory, number of components, accuracy, etc.), 
where a choice exists, he most convenient methods are those 
with which the user feels most familiar. It seems, however, 
mare sufficiently large analog computers do not exist, to per- 
mit study of big or strongly nonlinear systems. In this case 
the only alternatives are the hybrid and the digital computers. 

Table I synthesizes the main characteristics of the dif- 


ferent forms of modeling just described. 


46 








-szoqandwoo boTeue [Tews 10 AzOWSeW eZTS [Tews sortnbosay 
*STsou STYyA FO sdods 9yuA JO Ano st AAOSYA SUL 


"uoOTRIeAUSWSTAWUT 3[TNO TF Ia 


*3ndyno FO sonTeA uTW puUe xew 9Yy3 SeAZOsSoId 
yOTUM uoTSUedxS AONpoOAd sATUTJFUT ouA st stTduexs uy 
*“SqUTOd Mey 3e PeATSSep ST UOTANTOS SYA UsUM TedOTIOeIg 


*oOwTI Zegnduods ut A3soo ybty eyA st yreqnerzap AoLCeu suzy 
°“s,°o°p’d xoetTduoos Azan OF sTqeot{Tddy ‘*aeReanosoe AAdA 


Azessooou sowooeq [7/1 °*d “gt *Fey] AeaATSSqo 
ue JO OSn 3YR ‘SSATRICATASP FSATF SUA ACF SUOTATPUOCD 
TeETATUT UMOUXUN YATM SsuOTAeNbS sAeM OF potTtdde FT 


* (UOTIeWAOJ suery 
UTYLOTeH-AOUQqN, oy} *°a°*T) AeTndod sezow ATHhutseexzdsut 


*pesn 010 SOOUSTOFFIP OWT PAeMIOF JT STqQey3sun 


*butxeTdtjartnu AoF 
poou e st o79yq AeSuUuTTUOU ATHuozy3Ss St wa qsAS 3UA FI 


S}uUdUUIOD 


Auy 


suoasds 
uOTSUSUTpP TU 


Auy 


Auwy 


oOTToqzedAy 
pue OTToOqereg 


OTToqzredAy 
pue MOTT 


OTToqeied 


Co ts Gad FO OAT 


PDOPUSUNOD oY 


Auy 
Auy 


bo Teuy 
Oy Dt tC AH 


Teatbtd 
tO) Dia agar 


pes ona 
Te3tbhtd 
ZO ptaqAy 
HboTeuy 
TO Plaga 


HhoT euy 
TO ame Ar 


boTeuy 


xoanduod 


S,°ad°a°d ONILVWIOWIS JO SGOHLEW LNEYXEAATIGC AHL AO NOSTYWdWOO 


I Wiava 





OTAWO 
- aL NOW 


LLSSL 


LLSO 


LLSd 


LdSd 


LaS.L 


LOS 


LdSO 


LOSd 


pouzew 


47 








ZIT. -OPTIMAL CONTROL OF UNIDIMENSIONAL DISTRIBUTED PARAMETER 
SVouE vs 
foe LNTRODUCTION 

In the sequel of the general lines formerly enunciated 
with this chapter it is intended to create a physical feeling 
for the optimal control of distributed parameter systems, by 
summarizing different types of approach to the problem; for 
Simplicity only unidimensional systems will be considered. 

To accomplish this objective the solutions of five of the most 
representative problems are explained in depth. 

As mentioned in Chapter I the present work will not be 
dealing with derivations of optimality conditions or with the 
Memmematics of sensitivity, identification, stability, control- 
lability and obServability, although making superficial ref- 
erences to them. Each one of these topics has in itself the 
potential of an almost unlimited number of Scientific 
dissertations. 

According to a supposedly universal acceptance in this 
chapter the control function will always be represented by u 
fi@metne Optimal control by u*. Similarly, all the optimal 
quantities (states and cost) will have a * as Superscript. 

i the Open-Loop Control 

This has been the way most of the problems were worked. 
It requires the knowledge of the initial conditions and one 
Specific field of application is in problems of known 


meajyectory. 
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Schematically, the system representation is as fol- 
lows, where the double lines represent distributed quantities, 
generally in vector form for one-dimensional systems and in 
matrix form for multidimensional ones. The optimal control is 


computed and stored in order to be used when necessary. 


 Blatrey al 
Knowledge 
Of the syste 


Gemini zauion 





Figure 3.1. Open-Loop Control 


Motice that now the control is a function of the states. 
With exception for the requirement on the Knowledge of the 
initial conditions, all the statements in the last paragraph 
also apply here. 

For linear systems with quadratic cost functions the opti- 
mal law is linear and an equation equivalent of the Ricatti 
equation for lumped parameter systems may be derived and then 
solved in closed form. 

The cost is a well-defined function which can be computed 
in several ways, namely by numerical methods. 

The only input to the system is the control action, which 
may be included in the boundary conditions or in the p.d.e. 
describing it. At this point any of the techniques described 


in Chapter II must be used. 
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The subject of the realization of the optimality 
conditions was left open. Not much of practical value has 
been said about such a topic; however, a large amount of 
theoretical developments involving Calculus of Variations 
and Functional Analysis is available in the literature. 

Once again the approach may be divided; consider 
direct and indirect methods. 

a. Direct Methods 

The Direct Methods do not make use of the neces- 
Sary conditions for optimality. They consider very often the 
change in the cost function (6J) resultant from a small varia- 
mon (6u) and, starting with a trial solution, they follow 
optimization techniques, such as the gradient method, until 
§3>0 for 6usc. 

If the p.d.e.'s describing the system are reduced 
to ordinary differential or difference equation by techniques 
Such as described in Chapter II (Space quantization, time and 
Space quantization, eigenfunction expansion, transfer function 
approximations, etc.), then the optimal control problem is 
reduced to a lumped parameter, one for which there are well 
known techniques. 

b. Indirect Methods 

The indirect methods may also use the known re- 
Sults of lumped parameters optimal control. However, either 
uSing these concepts, or deriving other results, they always 
require the knowledge of the necessary conditions for opti- 
mality. This will originate equations such as the Hamilton- 


Jacobi~Bellman, principles such as Pontriagyn's maximum 
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principle and operational procedures such as the moment meth- 


oa® and dynamic programming. 


: Some other types of approach are possible, as for 
example the one shown by Balakrishan [Ref. 6] which includes 
constraints in -the form of a p.d.e. in the cost function. 


2. The Closed-Loop Control 


tts general Configuration 1S sheen etig. 32. 





Optimization 
Procedure 


ERGUKe 93.2. sheeadback Contre. 


Be. PHYSICAL EXAMPLES OF OPTIMAL CONTROL OF D.P.S.'s 
1. Heating Problems 
Suppose a large heating furnace used to heat metal 


billets is given and that these billets move uniformly through 


° This method was formulated in general terms by Kreyn 
in 1938 and used by Butkowskiy. It is a good method within 
its scope of application: linear systems with known 
eigenfunctions. 
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the processing zone. The requirement imposed on the plant, 
besides the cost, 1s that the outlet material must be at con- 
Beant temperature. In order to be able to Himpllement a good 
control system the engineer has to take in account the time a 
given section of the material takes to cross the furnace, the 
temperature distribution, the physical properties of the 
billets, as well as its thickness, etc. Intuitively, the pos- 
Sible control laws would take the form of changing the speed, 
or the temperature, or both. 

There are numerous constraints to be taken in considera- 
tion, as for instance in the furnace temperature T, (x,t), such 
that C,<Ty Cx ies Cys where Cy and Cc. are constants and x is 
the distance from origin. Also the interaction between adja- 
cent zones of the furnace must obey the inequality 


dT, (x,t) 
— $C , the billets temperature T. (x,y,t) must be less 


than a constant (here y 1S meaSured in the direction of the 


| — eae) 
thickness) and, Similarly, <a A : 

2. Heat and Mass Transfer Problems 

One example of this type of problems is the process of 

Seveng Of a moist material in a continuous drying unit. In 
this case the control action consists in compensating the 
changes in the composition of the moisture, porosity of the 
Material, layer thickness, velocity of the material inside 
the heater, etc. The constraints may assume different as- 


pects, a typical example being m<u(x,t) <M, [22 <c, and the 


temperature inside the material also less than a constant. 
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3. Separation of Liquids and Gaseous Mixtures 


When using big reactors their distributed nature must 
be considered and once again the D.P.S.'s techniques are re- 
quired. The problem may be formulated in three different 
ways, in terms of the control action: 

(i) Control of the source heat to obtain the highest con- 
centration of the substance, assSumimg a constant output 
atte; 

er) Control of the source heat to obtaim a maximum product 
output for a given concentration of the mixture; and 
(111) Minimization of the control costs for a given concen- 
tration and output rate of the product. 
4. Automatic Control of Large Hydroturbines 

In this case, the distributed parameter nature of the 
pipeline supplying water to the turbine has to be considered. 
The control of the turbine rotor is realized with a valve, 
but while this action should be fast, in order to minimize 
abrupt changes in velocity when the load is removed, it has 
to be limited by the resulting changes in pressure on the pipe- 
line walls, which cannot exceed the safety value. 

5. Gas Transfer Through Long Pipelines 

The problem consists in locating and controlling com- 
pressors along the pipeline such that the variations in the 
output pressure are Minimized. Considering only one compres- 
Sor situated at the origin (Single pointwise control) the 
boundary conditions will be of the type 


a (te) 


lI 


Dione) 


Viet) 


v(t) ' 


aS 








where p(o,t) is the pressure at x=o and v(2,t) is the gas 
velocity at x=%, This velocity equals vel Wise. lee pens 
a pmeeion of the load. 
The functional to be extremized is of the type 
Ah 
g= f |p*(t) - p(z,t)*lat , y21 
O 
The control of systems such as those just described 
is by itself a challenging task due to the complexity of the 
equations involved. Unfortunately, this task is made even 
more difficult by the fact that some of the parameters of the 
partial differential equation are not accessible, which im- 
plies the need for the use of parameter identification 


techniques. 


eee rLVE REPRESENTATIVE EXAMPLES 

A brief analysis of the contents of this section is given: 
Example No. 1 [Ref. 5] 1S an open-loop indirect method appli- 
cable to a large variety of D.P.S.‘'s for which a pointwise 
control is desired. It 1S based on a Space discretization 
technique that leads to a set of ordinary differential equa- 
tions to which the methods of optimal control of lumped para- 


meter systems can be applied. 


Example No. 2 [Ref. 5] has similar characteristics to the pre- 


ceding example but is a multiple pointwise control problem. 


Example No. 3 [Ref. 56] is the last of the open-loop problems 
presented here. It is also an indirect method and may be used 


Whenever the partial differential equation may be reduced to 
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the form Q(x,s)=G(x,s)U(s), where Q(x,s) and U(s) are respec- 


mavyely the transforms of the output and o£ Eke anpoue. 


Example No. 4 [Ref. 29] is one of the few closed-loop methods 
available in the literature. It applies to parabolic differen- 
tial equations with quadratic cost-function. The control is 


multiple pointwise. 


mempee NO. 5 (Ref. 28] presents again a closed-loop technique. 
It applies to partial differential equations with known anal- 
ytical solution and when the functional cost is quadratic. In 
the example shown there is only a single control but a multiple 
pointwise control is also possible. 


im Open-Loop Minimum Time Optimal Control of a Heat 
Transfer System USing Space Discretization and 


Pontryagin's Maximum Principle 


Consider the heat equation 


Z 


z Kod hee, sacl Ober ee) 
x 


= 


The boundary conditions are 


= eel _ _ 
Aoeleeg = OWE) = einen (3.2) 


where a and b are given constants, and 
= 0 (323) 


mime initial condition is 


q(x,0) = Aric. (*) (3.4) 


The pointwise control is bounded: |u(t) | <M 


S12, 








One way of reformulating a minimum time problem is 


to choose a convenient value for the final (t such that at 


¢) 
this instant the deviation from the optimal trajectory is 
Minimum. An adequate cost function is given by 

L 


J = f |q*(x) - q(x,t)| "ax , y21 (3.5) 
8) 


Discretizing equations 3.1 through 3.4, the following 


expressions result for L=N2: 


oy) Ee ia ie 


oe. Nema R72 
dt g 





dq. a 
eee iy R _ a 
a = k >, i =2,3,.--N-1 (3 La) 
A Gnt+i Gn | Fy IN-1 
Ange, SEW g 
elite g 
417I 
= en awe = b(u-q,) (3.24) 

q ~ q 

N-1 N _ 

=o 0 (3.3a) 


Meese is clearly illustrated in Fig. 3.3. 
The above discrete equations can be written in dimension- 


less form: 


> 


ie 


=e 2 — P& 2b 
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Gi, 
Co 
ae 
ae ee cr Fh at 
U 
ahs eee! (N-2)1 (N=) IL 
Bigure 3.3. Method of the Straight Lines 

dq 

1 Zi 
——ee — T — —— 
Te 5 ey hee Oe, £2 
O07 ; 
at = eal = 2g; + Sle 34) 1 P27 ce PN (3 lilo) 
Py 
qt ~ In-1 7° dy 

Peestituting (3.2b) in (3.16) the result is 

dq 

eA Me PO ee aes 
en oe ele 
dq; 
qt = ere: = 2g; + eel 1 Oe eee ele (3 ke} 
dq 


elen dy-1 7 Ry 


which is easily implemented on an analog computer. 


of 








The given p.d.e. may also be simulated easily 


remembering that 

i 
Ee rE (3.6) 
gives the voltage distribution in a non-inductive transmission 


line that has r and c as resistance and capacitance per unit 


length. Therefore the total resistance and capacitance are 


R=rL and C=cL. 





Figure 3.4. An Analogic Simulation of a p.d.e. 


The equivalent m sections circuit 1S represented in 


Fig. 3.4 where $= 28N , RC = 2 con 
1 
and the load R, = 


Return gebO.HGs 3.lc, ssuppose: Enat the following 


o because of Eq. 3.3. 


initial conditions are given 


Ay (0) 


(3.58) 


oa) | 4a Oeeeeet 





q, (0) = Tg (0) ( 


Equation 3.5 reformulated for the discrete case 
becomes 


eee. 
A= 1 


miach is to be maximized. 


(329) 


oe 








The system (3.1lc) is representable in the linear form 
as 
Sey ee EMGe ese) ae Io (Selo) 
and, aS a consequence, the necessary condition for optimality 
from Pontryagin's principle becomes also the sufficient 
condition. / 


Writing the Hamiltonian® 


28 2+38 


red eS Py (oyg U- FR) + do 


all. 
Le : : = : : - ° 
2,4 (ia 20549541) *Py (ya dy) Sa 
where the p's are the Lagrange multipliers. 
As it is easily seen, given that |[u|]<M , the Hamil- 
tonian has a maximum for 
u(t) = M sign p, (t) C3229) 


Looking now at A of Eq. 3.10: 





-5f 1 0 0 ----0 0 
1 -2 1 O----0 0 
DR ities ntat (3.13) 
0 0 0 0 —2 1 
0 0 0 0 Lt 


it should be noticed that it is symmetric and therefore its 
e€igenvalues are real. It also can be shown that the eigen- 


values are all negative. This last result shows that the 


U Reference 65, p. 234. 


: Reference 37, p. 188. 
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system moves to the origin as t+» and this implies the exis- 
tence of a unique optimal control. Because the eigenvalues 
are real and negative the optimal control will change its 
value n times, as illustrated in Fig. 3.5 for n=4. The ini- 
meee and final conditions are q, (0)=-1 and q, (T)=0 for 


1=1,2,3,4, 





Figure 3.5. The Optimal Control Law and the Correspondent 
OMhererbae 


2. A Simple Example with Open-Loop Multiple Pointwise 
Cone rod 


Consider the heat transfer equation 


og cs = 
ae , tL) vs + a(x,t) v(t) Te =P G) Wie ts) G3. leas 
between a medium at temperature u and a body moving with veloc- 
ity v20. The heat takes place between the x origin and x=L. 


The remaining information iS: a and v are known and 


q(x, 0) -= tic. (*) (SEs 


a0 Ae) 0 , Ocect (os 
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The problem consists in finding the control u(x,t) 
each that the functional 


at 


T 
3 = f |q*(t) - q(t,t) | ‘at, yet yey) 
0 


is minimized or equivalently, such that the time integral of 
the temperature deviation of the output material from a given 


law 1S minimized. The following constraints are imposed: 


(Ske) 


Discretizing the system inn linear elements, Eq. 


3.14 becomes 


aq. Gin 
1 aaa _ 
: ae a. (t) v(t) oa ao u, (t) (3.19) 





mie by Eq. 3.16, q(t) = Q and therefore it follows 





dq; u, 
Af a Bd. 4 2 oe te Big 1 A= eZ tess Soy (32.20) 
where 
= WE) See se) 
B = q and a, = a, (€) 7 (355.2 28) 


imorder to have the term with us fuec ot oOtnet eLmMe Lunes lons 


define a = a,q. and rewrite Eq. 3.21 as 
| 
| = ! a a (2222) 
oe eee 


mee adiscrete form of the constraints becomes 


C. S u; (t) < Cy 
(3223) 
Cy BS Sele ek. s c 


it] a 5 
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ama the cost takes the torm 


T 
g = f l[q'* - eee) acer y>1 (3E24) 
0 


The corresponding Hamiltonian is 


N 
H(p,q,u) = ye eee aCe) ert (Ray) |e) p, (Ba}_j+o,q!+u,) (3.25) 


i=l 


and, aS usual, the costate equations are obtained by 








dp (t) 7 : aE 
Gite 2 aq! 
or 
* 
m0 = & 
dt 
dp * 
pe CUCU BPs 44 pee = 12 ON (3.26) 
ai * 
n 1* eal : ' 1 
a 0 ( YPgia’* = elt! sign(q'*-qy) - Dy 
As it may be immediately seen, for Pp; = const. H 


takes the maximum value for a corresponding maximum of 


N 
F(p, ,u;) = ) p 


A us Satistying Ene constraints Inekq.s3.2 5. 


If L=22 it is possible to draw one elucidative two- 
dimensional picture showing the effect of the constraints on 


the optimal control. Equations 3.26 together with 


* 
a = ge mcie, cOnemol constrainus ance Hi(q=(0) 7 su (hiro. lig, 
T)6T = 0 2 lead to the complete solution of the problem. 


- Reference 37, p. 233. 
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Figure 3.6. Geometrical Interpretation of the Control 
Constraints. 

Butkovskiy [Ref. 5] describes two other methods used 
feeetnd the switching times of the control in minimum time 
problems. They are the method of harmonics and the method of 
parabolic approximations. Both give better accuracy than the 
straight lines approximation illustrated in Fig. 3.3 and they 
are mentioned here as a reference to the interested reader. 

Finally, since there are few references to multi- 
dimensional systems, it seems useful to refer to the above 
mentioned work for an explanation of a simple way of reducing 
the problem of optimal heat of a sphere or of a cylinder to 


the one-dimensional problem using the method of harmonics. 
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3. Sakawa's Solution of One Open-Loop Optimal Heat 
Transfer Problem 


Sakawa'ts solution [Ref. 56] is one of the most com- 
prehensive and complete papers in the field of practical ap- 
plications of distributed parameter systems. 


Given the normalized parabolic equation 


Z 
3 Q(x, =) _ og (x,t) 


for O0<xsl and O<t<T with initial and boundary conditions 


q(x,0) = 0 
a = afq(0,t) - v(t)} (3.28) 
Ble esA ie) a 

Ox bel = 0 


Mmrewecontrol u(t), O<t<T, is a fuel flow which controls v(t), 
the temperature of the medium. 
Equation 


y eee) Sapte) 2 ue (3.29) 


expresses the relation between u(t) and v(t) and is physically 
equivalent to a first-order lag. 
The objective 1S to minimize the functional 
1 


S(u(t)) = f {q*(x) - q(x,T) }*ax (3.30) 
: 0 


where q*(x) represents the desired temperature distribution 
at t=T. 

Applying Laplace transform techniques to Eqs. 3.27 
and 3.29 and solving for the given boundary conditions, 


mies) may be written as Q(x,S)=G(x,s)U(s). 
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From this equation, by the convolution theorem 


fe 
{ ofl ree bye at) cola: 
0 


ec) 


(3.2340) 
ic 
{ Gtk, c= 0) Win acm 
0 


Meee: All the procedures that follow may be applied not only 
to the heat equations, but to any partial differential equa- 


mfon reducible to the form of Eq. 3.3l. 


Rewriting G(x,s) as the Gquotrent ~Gix;s) — ee 
ME 1S easy to use the theorem of the residues and obtain 


NSH Ss. } s.t 
ae (3ues2) 


cea) =e 
1=0 


i 1 
t 
M (s;) 
where the M'(s,) 1S} elev; lSinLwetsatwis Cle IIE) eee S=S. 
By the use of standard minimization techniques des- 


cribed in detail in the reference paper, from Eq. 3.30, the 


optimality condition can be derived: 


il 
{ {q*(x) - q(x,t) }g(x,T-t) dx = 0 (33235) 
0 
1 
Metining f(t) = {| q*(x)g(x,T-t) dx (3.24) 
0 


memcomes out, by bottom Eqs. 3.31 and 3.33: 


dk ile 


er) a eu lj ey a (Xk, 1-1 )g Gz, fn )idxdy (0<1<T) (3235) 
0 0 


where up is a time dummy variable. 


The integral on the right 


i 
y(t) = f g(x,T-t)g(x,T-1) dxdu (3. 36) 
0 
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is a symmetric kerne1?? and as such has some important 
properties: 
— . . lal 
(i) There exists at least one eigenvalue, he #0 3 
a 
>, (T) = hs f y(T,H) >, (UY) du where 0..:(7) (0<1<1) 16) the 
0 = aa 
eigenfunction corresponding to ds 
(21) The eigenfunctions are mutually orthogonal, i.e., 
at 
J b(t) o;(t)at = 0, ford; #2, 


(111) The necessary and sufficient conditions for existence 


om SOlUELoOnS Of BQ. 3.35 12S that ) Ae c* must be 2con. 
1=1 
Ae 


vergent, where C. = f £(t) >, (t) dat ~) The optima! 
0 


control will be given by: 


oo (ie) = A.C. o,(t) . (3.229) 
al 


From Eq. 3.37 it can be seen that in order to obtain the 
optimal control it iS necessary to solve also Eq. 3.35 which, 
in general, is a hard task. Also it happens that sometimes 
Beemcolution of Eq. 3.35 does not account for control restric- 
tions. For these reasons Sakawa develops a numerical integra- 


{ 
Mme Of Eqs. 3.31 and 3.32 obtaining 


Z 


4 Z 
Sb) ORC Sees ier ibis, (3.38) 
i=o 7 1 3j=0 ae, 
10 See Ret. bl, so. 25, Lor detailed definition sana 
properties. 
11 


3 means “Such that." 
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where the Cc; 's and ig are obtained in a straightforward 
manner. <A way of minimizing Eq. 3.38 is by quadratic program- 
ming, for which Sakawa gives appropriate references. However, 
because a linear programming technique is much easier to 


implement, he considers another cost function: 
1 . 
WG ea) emeeiere (<2) | cts (3.39) 
0 


From this he derives, in a similar fashion as for Eq. 3.38, 


the following equation: 


1@l * N 
Tue |) Cela. = o 21, Ul | ; (3.40) 


The control constraints are taken in account in the linear 
fer nonlinear) programming. 
In the sequence of this last procedure the author ob- 


tained the following curves for the parameters 


o= 10, y= 0.04, g*(x) = .2, #4x-n-= 20 

Summarizing, Sakawa's method is considered to be quite 
Valuable for single pointwise control of unidimensional sys- 
tems. Although it iS not impossible to generalize this meth- 
od for multiple control using multiple Laplace transforms and 
the superposition principle, this seems to be a huge task and 
therefore the jachod 1S not recommended for such cases. 

Quite recently, Chang [Ref. 9, August 1970] developed 
one algorithm based on a modified steepest descent method and 


solves the same problem as above obtaining very good agreement 


With only two iterations. 
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Meeuce 3.7. Optimal Control 
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Mere 3.9. Optimal Control 
(T = .4). 
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Figure 3.8. I. Temperature 
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II. Temperature Distribution 
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Figure 3.10. 1. Temperature 
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4. Greenberg's Solution of a Closed-Loop Multiple 
Pointwise Optimal Control, Problem wi theo madu ac ie 
Cost  Punceren 
a. Pointwise Optimal Control 

The mathematical derivations used by Greenberg in 
his paper [Ref.29] are hard to follow if the reader does not 
have a deep understanding of Functional Analysis and theory 
of operators. However, the consequences of this method are so 
meeortant and its implementation 1S so simple; if compared 
with the complexity of the derivations, that it is worthwhile 
to analyze its contents in some detail. 

The typical approach to the problem represents a 
distributed control law by a finite number of lumped controls, 
but this is not applicable to a feedback control. Also the 
modal control, by the reasons already pointed out in Chapter 
I, 1s difficult to realize if the eigenvalues are hard to ob- 
tain or if they converge slowly. As it will be shown, the 
present method does not suffer such drawbacks. 

The author starts with the knowledge previously 
demonstrated that systems described by parabolic differential 
equations with quadratic cost functions have an optimal feed- 


meen. control. 


Writing the parabolic differential equation in 


mae form: 
me) = Aq(t) + B(t)u(t) , alo) = qos (3.41) 
where A is a partial differential operator (generally of the 
42 
type 77) and B(t) is a time variant control operator defined 
Xx 


mee) < t < T, the functional cost is defined as 
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st 12 
Ab 
S= {| [<Ole)aq(e) (ate) -+u- (2) RiGee (3.42) 
0 ~s ~s ad 
R{t) 1S a kxk positive definite matrix and Q(t) a positive 
@eperator; k represents Ehe number of desimedvcon ene Mie omuecnr 


After non-trivial manipulations the optimal con- 


m@eol law is found: 
ut(t) = - R7(t)B(t)f k(t,c)q(t,z) az (3. 43) 
_ E < D 


where Be Gepsy = K(t,X,x,), 1 = 1,2,---k with K(t,x,z) the solu- 

tion of a given integro - differential equation, and 

B(t) a kxk diagonal matrix, with easily computable elements. 
The block diagram correspondent to the above men- 


tioned description is represented in Fig. 3.11. 


, K VeCeorm 


u(t,x)=B(t)u’ Dynamical System} q(t) 


nye ; GqGiSsttribured 
ea Aq+u(tx) 


quantity 


‘NL 













Pointwise 
Control 
Operator 





ROU) | k(t x XJdx 





Bmeoure 3.11. Block Diagram for Greenberg's Method 


y Lie Symoolnw <x," Meals Inner Produce Of the vectors 
Paeend y. 


ne 
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B(t) is defined using spectral decomposition of operators, 
a subject that will be treated in detail in the next chapter. 
b. Infinite Terminal Time Problem 
The above description can be quite simplified in 
the case of a terminal time problem because the operators Bo 
o, R and K the feedback operator become time invariant. 


Re Me sent 


ll 
OC dna) VO) (3.44) 
Boers V(x) is the vector of the first N eigenfunction of A 
and Q aN x N positive definite matrix. Define K as the posi- 


mye definite solution of the matrix Riccati equation 


Mme - KA + K Vo - 


BR = Baw kK oO. =— Ova A aN xN diag- 


onal matrix such that Ass equals hs of A. Also define V as 
a k xX N matrix whose eee element is Vis = Wa asa Ic The optimal 


control is given by 


ES) 


i 
ve) 
ee) 


BV K fv(t)q(t,2) az 
: (3.45) 


where q(t) is the N vector of modal coefficients of the state 
distribution q (t,x) 3 

Due to Ba. 3.44, a2£ the Gpeimal icontrol existsere 
Seaomot be improved by feeding back more than N modes. An- 
emer important conclusion, is that whenever (On 1 On) is a 
positive operator it follows that Kye Ry 20 This can be 
Stated as "monotone approximation of the state weighting 
Operator originates monotone approximation of the optimal feed- 


mack Operator." 
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Greenberg finishes his paper with a very simple and 
illustrative example. 


at, 


5. Graham's Solution of a Close-Loop Pointwise Optimal 
ConeEro! Preplemn 


J. W. Graham's paper [Ref. 28] presents two different 
approaches to the solution of an optimal control problem char- 
acterized by a given parabolic differential equation with a 
frequent type of initial and boundary conditions. The author 
uses two types of approach to the problem; the first one deals 
with the matrix Ricatti's equation and the second one con- 
Siders Kalman's equation, together with root-locus techniques, 
SPemstituting, in a certain way, like a smooth transition be- 
tween the classical and the optimal control problems. 

Some details in the derivations are very hard to fol- 
low due to gaps in the explanatory theory; for this reason 
additional information will be given in an attempt to clarify 
the points where more omissions were found. The nomenclature 
used is consistent with the remaining of this chapter, but the 


reverse of the one in the paper regarding the control and the 


temperature. 
The p.d.e. 
: » BG 
SI (x,t) =v ed tc 05a; O0<xX<a (3.46) 
ou 0,2 


Morgiven with initial conditions 
q(x,0) = Q. (x) (3.47) 


eee boundary conditions (B.C.'‘'s) 


dg 7 
ae (Opel = 0 


Ked(a,t) = £(t) - q(a,t) = u(t) (3.48) 


WZ 








where K is the thermal conductivity and f(t) is an external 
available temperature input. All the other symbols are de- 
fined as before. 


The Ycose fUncGtelon .uSae fincas 


© @) 


3 = f [q(a,t)-q.,(a)]* + Rfu(t)-u,.)7at (3.49) 


; S 
where the subscript ss represents steady-state and Ris a 
positive constant. The physical meaning of the equation is 
obvious. 

phe analy ele SselurioOn Or the pide) for etiesea wom 
Beundary conditions and piecewise constant control (u(kt)) is 


given (although not necessary) : 











2 Zone oe n Deeg 
eax , t) te Wee S2gaeg — Ay (-1) ier eee n1x 
ak) Ree eS L 5 exp | 55 7 COS Same 
6a tT n=ln a 
(3250) 
meee tkKT<t<(k+1)T]. 
After this the author states: "The original sampled- 


data equations in discrete form are transformed to the normal 
Or diagonal form in a straightforward manner and from these 
equations a continuous-time set of equations is readily ob- 
tained. Hence the system equations may be written as the 
following: 


W(t) 


PINE te) Sees UT (ie) (325.9) 


O(E) = C Wit) | (r529) 


where W(t) is an N dimensional column vector, A a NXN diag- 
Onal matrix, B and N-dimensional column vector, Q(t) the 
temperature of the system at x=a and C is an N-dimensional 


mector."™ 


de 








The way Eqs. 3.51 and 3.52 were obtained does not 
seem So obvious as the author pretends. One possible tech- 
nique that can be used to find them is the Bubnov-Galerkin 
method described in Appendix A. 

Back to the original problem, the ftunctrvonaley sean oe 
Earther simplified by defining 


Q*(t) = Q(t) - Q.. 


(3.5.3) 
* = ca 
u* (t) Gt) ees 
from which 
W(t) = AW*(t) + B u(t) 
_ oy (3254) 
OG) =e ws (te) 
meth initial conditions 
O*4o) = =O 
ad (255) 
We (Oo) sia 
SS 
and it comes out 
~ 2 2 
J = f Ors a tee ect (3.56) 
0 
or 
a T 2 
J =f (W*(t) SW*(t) + Ru*(t)“]dt (357) 
0 7 aoe 
where 
C2 Che c 
1 ie) Salon 
Cac on c 
ee | ea 2 --- ©2% (3.58) 
' ' ‘9 
Co Cua 7 Cu 
Semecne Ci's (i=1,2,----- N) are the elements of the C vector. 


Note that in the above equations the superscript (*) is a 


mean for identification of some variables after 
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transformations on them, rather than being the usual con- 
wemtiton fOr optimal function. 
: At this point the author considers two different ap- 
Beoaches to the computation of the optimal control 
a. Matrix Riccati Equation Method of Solution 
pataune of the linearity of the given system and 
also because a quadratic cost function was chosen, it is 


fm@ewn that the computation for the optimal control may be re- 


mreea to the solution for Dol the arceaci seq Watt On: osama Fr 


mee? 52-7/5) 
dl 
: P(t)B B P(t) 
P(t) + S - as as = + P(t) A + AP(tE) = O (3. 5.00) 
meen boundary condition P(t,) = 0 peli s beundary yeond1t v6EM 
Pemes from the definition of J. If J had another quadratic 


term outside the integral, i.e. the term wei owe, the boundary 


condition would be P(te) 


The optimal control law is given by 


= G. 


ob = 


* => —_ 
u (ay ee R B PW(t) (3.60) 
where 
ey PN Sey O61) 
- eee © 


whenever the system is controllable and has P(t, = 0). 
P is obtained as the solution of 


_ a) 2 Nes Ban uae = 6 2 oO, (3.62) 


Peeby SOlving Eq. 3.59 backwards in time from the known 


condition P(~) = 0 until convergence is reached. This last 


12 








method can be implemented very fast in a digital computer 


ioing the algorithm 


P(t+A)=P(t)+A{-P(t) A - AP(t)+P(t) BR /B P(t) -S} (euce 
The graphical interpretation of the procedure 


iemogiven in Fig. 3.12 





O at = i 
ST SSE State transient 


fMmeedre 3.12. Graphical Interpretation of the Solution of 
Re uaevenes. 6 7. 
It 1S important to note that P(t) is independent 
of the states and can be computed once J is specified, before 
the optimal system is started. 


Once u*(t) is obtained from Eq. 3.60, the final 


pe 
result DME ys 1S given by BE) oat + — 
b. Kalman's Equation Method of Solution 

The foregoing method could be applied without the 
need for reducing the state equations to the diagonal repre- 
sentation. This is not the case in the present paragraph, 
because the nice results derived from Kalman's equation, and 
explained in deep detail in Chapters VII and VIII of Ref. 60, 


are much more easily used for the diagonal form. 
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Figure 3.13. The Global System 
Suppose an arbitrary i row of ann dimensional 
state equation: 
w. =A.W. + b.u. (3.64) 
1 a Mego Y a 


Taking the Laplace transform and ignoring the 
iMrtial conditions for the sake of clarity, the distributed 


system can be represented as 











_ Di 
S-)y 2 
b 
SS) Cy 
Us) ! (+ ais) 
| ne. ie 
1 S- L 
| es ean See | eS 
Da 
S=)y 
Begure 3.14. Approximate D.P.S. Transfer Function 


(becomes exact as n>) 


77 








N-1 N=2 
Q*(s) _ a> + AS + —-—--— + a5 


G(s) = U*(s) = =(S=\nGS] |. ==aitc> a (S-h,)--- (8) (3256 5) 





But, as shown in Ref. 60) £eem Ea. 3.66 5alse 


known as Kalman's equation it can be derived that 


{1 + G(s)H, (8) |* = panes 2 (3.66) 
where 
l1+G(s)H__(s)|7 = (14+G(s)H__(s)) (1+G(-s)H__(-s)) 
eq eq eq 
2 (30a) 
[G(s) |“ = G(s) G(-s) 
and that the optimal control is givey by 
opt = r-hW*(s) (r=0 in the present case) (3.68) 
such that 
* * See * 
a.) - hW* (s) hywi(s) + h,W3 (s) + + Aywy Cs) 
eg aoe * 
QO* (s) cy wi (s) + cow (s) + + Cywey (S) 
(Sim69) 
This is easily illustrated in Figs. 3.15a and 
Bei5sb. 





Figure 3.15a. Block Diagram of the Compensated system 
(r=0). 
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-_s 
Y) ! 
am 


aa 
aS 
nfl [a] fre 


— =e =e se ee ee eee ee 


= 


Figure 3.15b. Another Configuration of Figure 3.15a 


In order to Obtain the eGllementcs of h, solve Eq. 


3.66 or equivalently, the equation?? 


Ete) sees 


InG(s)H. (5) = [1 + R ] (3.709) 


where the [ 12 ({ ] ) indicates the poles and zeros of 


G(s) G(-s) 


it R 


Ehatedae inthe eft (rregit) siakt so planes 
The denominators of both sides of Eq. 3.70 must 
be equal and the numerator of the left-hand side can be ex- 
pressed in terms of the h;'s. Therefore, 1t is only neces- 
G(s) G(-s)) 
R 


wemy to find the L.H.P. roots of [l + either by 





13 


Fe G(s)G(-s).* 


G(s) : 
1s 


V R 
From Eqs. 3.70 and 3.71 it is possible to write 


G(s)G(-s))_ 


ies R 


Plt 








G(s)G(-s),~ G(s) G(-s) ," 
R 


ep) (Ss) 1 litG(-s)n (os) -[1t }) [1+ R 


Assuming that G (S) has poles only on the L.H.P. it can be 
Shown that the poles of RES) 5S) are also on the L.H.P. 


G(s)G(-s))" 


and therefore TANG WSL VS) = {1 + R 


a) 








root-locus or factorization techniques, after which the coef- 

ficients of equal powers of s are equated. 

The results obtained by Graham in two given prob- 
lems, by this method and Riccatti's equation method, are in 
very good agreement. 

Pee SD ELORT COMMENTS ON OBSERVABILITY, STABILITY, CONTROLLABILITY, 
NONLINEAR PROBLEMS, STOCHASTIC CONTROL AND PARAMETER 
POENTLIFLCAT ION 
The purpose of this paragraph is to give some references 

to the topics enumerated in the title. 

Because Ref. 52 furnishes all the necessary information 
about the references prior to the end of 1969, only those 
published after this date will be considered. 

1. About Observability and Controllability 

Not much has been published in this field, and so far 
as it is known, all the reports dealing with the observability 
problem are already mentioned and superficially analyzed nue! 
Bei. O2. 

Also the controllability problems have been in a peri- 
od of abandon until recently (1970) when the International 
Journal of Control published a paper by Herget [30], written 
in 1968. The author considered a distributed parameter con- 
trol problem of a linear system with constrained control, 
either distributed or at the boundary. 

The system studied is represented by the equation 


dQ (x,t) 


aE = A(x) Q (x,t) + £(x) u(t) ; XEQ ce L CST) 
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with U(X) Q(x,t) = Q for x in the boundary, and Q(x,0) = 0 
LOX xell. 

In the above equations A (x) Zepeesenivs jap areraiicne 
ferential operator defined relatively to %, U(x) for boundary 
x 1S an operator denoting the boundary conditions and E Coa s 
a function also defined in 2. 

The paper shows what are the reachable states of the 
above type of system for two classes of problems, namely with 
and without magnitude constraint in the norm of the control 
( ul] 7). 

wee The Stability Problem 

In the field of stability, some recent papers seem to 
be of great interest because of the wide variety of systems to 
which they apply. 

The following references are mentioned here: 

(1) Kathri [Ref. 35] uses multiple Laplace transform tech- 
niques in partial differential equations which can be written 
in the form 


i yi ma yk+5 M 
=0 *9tt k=l j=l *'J 9t*%dxJ me 





gy 
1 ne ox 


= U(X, ec) (3272) 


= 


for a large class of boundary and initial conditions. 
The system was reduced to a transfer function expres- 
Sion to which a simple stability criteria was applied. This 
report is quite comprehensive in the sense that the author 
Birnishes a step by step explanation of all the details. 
Later, the same author published a new paper [Ref. 36] 


in which he considered the same type of systems as before but 
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applied to sampled data problems where a memory-less, non- 


linear feedback was present. 


mer) Ansari [Ref. 2] studied the sufficient condition for 


Stability of a system described by 


0Q (x,t) SOS) gq >. = 
me to U(t)——— 2 = — Bu) 2 (x, £) aoe ee} 


Bob O<x<L , H>O and u>0 , where the control u = £(T(LZL,t)) was 


Siem that aS and O(0,6) £(0(L,t)) are <ontanueus anedenone— 





tonically decreasing functions of u and T(L,t) respectively. 
This condition can be verified very easily and although ap- 
plicable only to a small class of p.d.e.'s, englobes some im- 
portant diffusions systems. 
(i111) D'Souza [Ref. 14] treated a much more general problem 
than those mentioned previously. He considered a class of 
causal dynamic systems represented by nonlinear vector-matrix 
p.d.e.'s. For a good understanding of the developments car- 
ried out on the given reference, a good background in Func- 
tional Analysis iS necessary, with special emphasis in 
Green's function type of problems. Because such theory is 
not yet within the range of the knowledge of the average engi- 
neer, it seems to be useful to rewrite the paper in more 
practical terms, if a larger audience is desired. 

(iv) Kastenberg [Ref. 37] derived the stability conditions 
for nonlinear feedback control systems described by parabolic 


p.d.e. of the following type 


mee oY 855 Sx, tL Oy teat CON PEO) 
ij + 1 °*5 i : i = 


(3.74) 
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where F(x,t,Q) represents the feedback control law, which 
must be negative definite. The initial and boundary condi-— 


=rons are 


Q(x,0) Qo (x) fOr xX imside Lets tue gone. ide tani ona 
Q (x,t) = 0 for x in the boundary region the type of 
system studied is cf wide application and the results derived 
seem to be easy to implement. 
3. Nonlinear Systems 

Although published in 1966, the paper by Uzgiris and 
D'Souza, "Optimal Control of Distributed Parameter Systems 
with Nonlinear Boundary Conditions" [Ref. 67] SHG not receive 
much attention. This reason, together with the fact that the 
procedure developed is relatively simple compared with the 
complexity of the problems it can solve, motivated the present 
geLagraph. 

The authors considered the case of a one-dimensional 
linear heat equation with nonlinear boundary conditions which 
were made to include the bounded control function. Then, by 
discretizing the space variable, a set of nonlinear differen- 
tial equations was obtained and reduced by Laplace transform 
and convolution techniques to a nonlinear vector integral equa- 
=mon, £rom which Ene optimal control could be derived. Follow- 
ing standard methods, the final problem was reduced to the 
solution of a Volterra type equation which can be solved by 
known procedures. 

4. Stochastic Problems 
This paragraph considers the recent papers by 


Tzafestas [Ref. 66] and Seinfeld, Gavalas and Hwang [Ref. 61]. 
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The first one deals with the identification of para- 
meters of D.P.S.'s by reducing the system to a lumped para- 
meter system for which there exists known parameter 
identification methods. The last paper contains the deriva- 
tion of a nonlinear filter for systems with disturbances in 
the initial and boundary conditions. 

5. Identification of Parameters 

The present chapter is concluded by mentioning the 
modern method developed by Fairman and Shen [Ref. 15] in the 
identification of D.P.S. This method, named by the authors 
as the "moment functional method" is applicable to one- 
dimensional heat or wave equations and in the case of the 
heat equation it may be extended to the problem in which one 


coefficient is a polynomial function of time. 
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IV. MODAL CONTROL THEORY 


The theory of modal control was first developed by 
Rosenbrock [Ref. 53] for lumped parameter systems and the 
concept generalized by Murray-Lasso [Refs. 25 and 45], and 
Seuld (Refs. 25 and 26] in order to include D.P.S.'s. 

The above references give the general theory of what is 
going to be done in the rest of this thesis. This theory, 
although applicable in many Situations, suffers from several 
drawbacks. The major ones are the need for computing the 
eigenvalues and eigenvectors of the system's operator, calcu- 
iirwom of the inverse of Singular Matrices, need for © con 
Peeling Clements 1£ r modes are to be changed, difficulty of 
application in the case of repeated eigenvalues and also 
limitations when the modes are complex conjugates, in which 
case only the real parts can be compensated. 

Some of the inconvenient parts just mentioned were 
avoided by Simon [Ref. 62] and Foster [Ref. 18], and they 
furnish very interesting reSearch material in the optimal 


eomecol field. 


A. GENERALITIES 

In Fig. 1.1 is shown how a pointwise input produces a 
distributed output. The idea of the modal control consists 
in uSing the meaSurements at different points, uncoupling 


these measurements and obtaining the coefficients Ws of the 
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eigenfunction expansion of the output. 74 Then, compensate 
these eigenfunction's coefficients, compare the compensated 
values with the reference Signal and build with this informa- 
tion the adequate control law. This is done in a synthesis 
procedure conceptually identical to the analysis procedure 
that gave the coefficients Ws. The global picture of the 
control system for a heating problem is indicated in Fig. 4.1 
and it shows a baSic assumption that was made: The bandlim- 


Meedness—~ of the three functions considered here, namely 


14 It is shown in Ref. 45 that in the case of a system 
described by an operator L such that the output is 
y(x,t) = Lm(x,t) and 


-] : 

a) L exists, 

b) The eigenfunctions of L are products of functions 
of time and distance, 

eo) § tS completely continuous in x, Eat aS PO Saye 
has a purely discrete spectrum (eigenvalues content) and the 
only possible cluster point is the origin, then it can be 
stated, using the concept of generalized Fourier series, that 


y(x,t) = } w, (t)u, (x) 
dL 


mA te) ) He (t)u, (x) 


Al. 


where u. (x) is the ad SLgeniunecereOnmeOor f,NObeMeceooramn, 
Sine or coSine functions, and forms a complete set. 


Also the operator L is such that its adjoint, L*, has a 
complete set of eigenfunctions Vv, (x), Wihken 2S OrEnOgonc tae 
tier x) . 

aL 


Using the orthogonality properties in the first of the 
above equations and assuming normalized distance and eigen- 
functions it follows that 

aL 
ws (t) = f y (x,t)v, (x) dx 
0 


15 By bandlimitedness it is meant that the functions are 
considered accurately described by N eigenfunctions, instead 
Of the infinite number that exactly characterizes all the 
partial differential equations. 
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reference temperature, output temperature and control 


function. 


B. LUMPED PARAMETER SYSTEMS 

The baSic idea in the modal control of lumped parameter 
systems is the independent shifting of the lower order eigen- 
values such that the speed of response of the system is 
increased. 


In the case of an uncontrolled linear system such as 


ixe 


=Ax (4.1) 
it 1s possible to obtain x as a linear combination of terms 
containing the eigenvectors and eigenvalues of matrix A. For 


~~ 


this the transformation 
y=H x (4.2) 
is made, from which it follows that 


Ht y (4.3) 


~ 


wh 


= H 


tp 


If the eigenvalues of A are distinct it is possible to 


eotain 
Ay 0 (es 0 
PeeHe =o 0 A, Ooo. 10 
| 0 0 24 dn | (4.4) 
and therefore 
het 
eae aoe (Aes) 
where a, = y; (0). 
inal ivieshirom EG, 4542 e)ReEt, 46, p. Liz) 
N Ast 
—— ) Otel: 16 (4.6) 
ee dis oa 
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| [-] woodx pyrey (x) 


(XS) “{ ) mx) x,t), | 
Analyzer) Oe pis! Syn thesiZer =a ee 2 eu ee 





Figure 4.1. Theoretical Implementation of the Modal Control 
ye el (DAIS is 


Here u, is the ee eigenvector of A and hs is negative if 
the system is Stable. 
The purpose of the modal control is to shift the eigen- 


values in such a way that finally the following relationship 


is obtained 
H Cl ese) 2 
ee cu ce aaa (4.7) 
= se dea. 
with kK; a negative real number and as the distribution of the 


initial state into’the ae mode. The way this result was ob- 


tained 1s explained next. 
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1. Ideal Case 
The system is desc 
x = 
and the output is 
ae 


~w 


in this, adeal= cace 


mroethe dimensions of A an 


ribed by 
Avxs+ Bau (4.8) 
C x (ao) 


the dimeneitons or sand wearer 


aB are NxXN. 


~~ 


From basic Linear Algebra it is known that in the 


case of square matrices with simple and real eigenvalues the 


eigenvectors of A are orth 


~~ 


ena both matrices have the 


From the above sta 


Au. 
eich 


Nw 


ik 


and representing now U = [ 
the result is 
AU= U 
VA= A 


Wien. ene (dit he rent 


eouebe written 


—_— 
=— 


VU it 


and from this equation and 


A= UA 


~e 


Writing now the st 


Mag. 4.3 


ogonal to the eigenvectors of Me 


~ 


Same eigenvalues. 


tement it follows 


AGU 
C4 105 
pr 
se 8 r 
Nal 
tl, ul ] and V=lv a 
il oe So ee |e 
S Ak 
L~N | 
A 
7 (4.11) 
Wie 
u,'s and v,'s ate normalized st 


(4 2} 
Ghe precedent Set it is Obtained 
ve (4543) 


ate equations as indicated in 
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Figure 4.2. Block Representation of an Uncontrolled Linear 
System. 


L =} 1 = ‘ 





K 


Figure 4.3. Block Representation of a Controlled Linear 
system. 
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ix e 


= CA Beye (4.14) 
Finally, because of Eq. 4.13, a convenient choice for 
B and C would be making them equal to U and V, which gives 
x= U (A+ K) Vv x (4.15) 
The resultant matrix has the same eigenvectors as A a 
and the eigenvalues have been increased as was desired. More- 


over, the shifting of the eigenvalues is realized without 


Mieeraction, due to the diagonal form of A and K. 


2. The Number of Manipulators (r) rd Is Less than the 
Order of the System (N) 


Having chosen C = V in Eq. 4.14 the same->equation may 


be written as 


Ix 


- YA+BR) Vx (4-16) 


Now, in order to have B, aNxr Matrix, again as a 


Square matrix, make 


<Y> <N-r- 
; 1 
K le r 
BK = [U_} 0] a ae G 
~~ os = tis ’ ati: a eae ae a 4 (4, 17) 
xr NX(N-4r 
ee veer : 7 


16 Au=rzX Uu 


ed 


Suppose the eigenvectors are different and call them w. 
Therefore, 


U(A + K) Vwe= (A + K) w 
Awt+Kwe= AwtkKw 
But, the eigenvectors of A are unique and so w = u 


ay 


A manipulator is the device to which the control law is 


applied. Its output is the physical input to the system (e.g., 
beat flow). 


onl 








Smt Law iy7 


4 

A 0 r 
U A= ( fUy gd | t (428) 

g | oe + 

and from these three last equations 
laser 0 
X= U ae +o imal er ee LS) 
\ 


which shows how it is possible to shift arbitrarily the r 
lowest eigenvalues. 

Reference 26, pp. 246-248, contains the theoretical 
development for the situation in which the eigenvectors 
Is 6 6s +o oe are not known accurately. 

The procedure consists in making 


B= D(V_D) ~~ = {s: 


Ze b_] (4.20) 


pPorece- = 

and the final results are similar to the ones obtained pre- 

viously, with the exception that disturbances in the lower r 
Modes cause also disturbances in the higher N-r modes. This 
is not always a serious drawback since the influence of the 

higher modes can generally be neglected. 


3. The Number of Sensors (m) is Less than the Order of 
the System (N) 


The derivation for this case is also in Ref. 26. It 
is a little long and does not follow exactly the same lines 
as the distributed parameter case. For this reason this 
paragraph will only be concerned with the basic idea of the 


process. 


a2 








Start forming an, Nxmamatwisc 


W = [u, vee UL ULL ee a by uae oon (4.25) 
where p is the total number of modes taken in account, 
fe—- M- p, and b is defined as in Eqe 4.20. 


Then, form a Nxm matrix E such that e, (i=1,...,p) 


approximates as closely as possible of u; (2=1,...,p) and 


on ~e 


ona b. — U; for the last q elements e; Ome 
Define F, nxm, such that 
Fe = (om) oa" (A227) 


This matrix may be written as 


F [f,--.f, Foon OSs Broo bee 


(4.23) 


[F iF , | FX] 
Cries CH ta 


and from it obtain C defined as 


rT 6 | «r-q> f+ 
Cc’ = [F_ + F*: F ]n (A324) 
~ ~q 2 | i-q cv 


which has similar properties to C as used in the ideal case. 
The major difference is that the rows of C are no longer 
orthogonal to the eigenvectors u; (i=ptl,...,N), corresponding 
to the neglected modes. Therefore, the control of the lower 
modes interacts with the higher ones not allowing a perfect 
uncoupling. Generally, under the assumption that the higher 
modes are sufficiently high, the interaction will be small. 

Figure 4.4 shows the block diagram illustrating the 
present situation. Next it will be explained how the fore- 
going concepts can be extended to distributed parameter 


a ocems. 
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Figure 4.4. Block Diagram when the Number of Sensors Is 
Less than the Order of the System 


See OLoTRIBUTED PARAMETER SYSTEMS 

Before starting reading this section one should become 
familiar with the contents of Appendix B. 

1. Computation of the Feedback Control Law 

Following lines similar to those that lead to Eq. B.9 
it results that in the case of m and y being both functions 
of x and t, the diagonalization procedure gives 

w.(s) = dA, (Ss) us (s) (4.25) 
a OAS) and Ba MB, are the Fourier coefficients of y(x,Ss) 
and m(x,s) and NBG, are the eigenvalues of the operator. 

The truncation of the higher order modes allows to 
describe the control system in matrix form as shown in Fig. 
aD 

Ca liaing W the closed-loop matrix and using the 


properties of matrix algebra and of uncoupled matrices it 


can be written: 
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megure 4°55. Block Diagram Of tne Closed—-Loopysyseem 


W= (LF + 0) ' 4 | (4.26) 
Peom this it follows 

WL = (LF + aes (Ay 2s 

LW = LF + I (4a) 
pee tinal ly 2 

wo = F + ie (4.29) 


The system can then be represented as in Fig. 4.6 and 


for each loop the result is 


Fin See z 
OS) = Ww (8) X(s) (4-30) 


where the o,(S)5 w, 's) and d's) are, respectively, the 
eeagenal elements of Ff, Wand L. By a suitable choice of the 
>, Ss 1t is therefore possible to make each closed loop be- 
have as fast as required. 

The situation illustrated in Fig. 4.6 although desir- 
able is not possible in practice because of the existence of 
sensors restricted to certain positions and of manipulators 


Bimmen Cannot give exactly the desired output distribution. ft 


will be shown next that in a physical system it is possible to 


2 











Fugure 4.6. 


The Uncoupled System 


obtain the coefficients w. from the output meaSurements, and 


a., the actual input to each manipulator, from the knowledge 


of the distribution produced by the manipulators. 


MEE DStermNination Of Ene Pourier CoecEELTeclents Of tme 
Output from the Output Measurements 


Under the assumption of bandlimitedness 


N 


a ae a u, (x) (Aisa) 


Similarly, the measurements at XprXore-++Xy May now be written 
as 
N 


y (x, ,t) e) w, (t) u; (x4) 
, 1=1 


: (A= 32) 
° N 
VCS 2) = a a uw, (X\) 
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Oia, IneMatai <<) comm 


Y= Ut (47.33) 
where 

y (x, ,t) Uy (x, ) uy (x,).. uy (x) w, (t) 
y =| Y(X5-t) , we U, (X54) us (x5). uy (x5) Q = | M5 (t) 
¥ (Xyrt) Wea espa ese oo othe al | Wy Ct) | 

It follows that Ss) 

ee 
eae te (4.35) 


and it 1s proved in Ref. 45, Appendix C, that there is always 
a set of points XprXgree Xp such that a exists. It was 
verified in the simulation procedure of the present work that 
even in the case of the sensors being relatively close, it was 
possible to reconstruct accurately the coefficients Ws. The 
only differences are in having higher values of the elements 
in ay and in the sensitivity to errors in the measurements. 
There are several criteria described in Refs. 18 and 45 for 
the minimization of the sensitivity. The one used consists 


in the minimization of the maximum eigenvalue of the product. 


( U )° U (4.36) 
and is valid whenever the number of sensors (S) 1s greater 
or equal than the order of the approximated system (N). 
The symbol (+) represents a pseudo-inverse matrix and it will 
be defined in paradraph d. If S=N, the pseudo-inverse 
coincides with the inverse. 
When S<N the sensors must be positioned according to 


the minimization of the norm x - i U (x, )| , where now U 
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is function only of one element of the set {x, }-u" is also 
a different matrix and it can be obtained as explained in 
Ref. 18, pp. 60-63. The minimization 1S equivalent Eo veb= 
taining the minimum value of the maximum eigenvalue of 

+ 


Pee OG 


3. Determination of the effect of the manipulators 





Call H the distribution caused by the ce manipulator 
and decompose it in an eigenfunction expansion. 


H, (x) = by Uy (x) thy guy (x) +... FD uy (x) 


: (As) 
Hy (x) = by ty (8) thy ots Ca) +. . - tb Uy CX) 


The instantaneous output of the ee Manipulator (assuming 
instantaneous response) for the actual input a, (t) is 


therefore 


a; (t) H, (x) (4.38) 


and given that the desired input 1s 


N 
cet Narre e) lie (50) (4.39) 


Equating both 


a N N N 
a OE ills amare ig 10 by cadet oo 
or 
Bia = 1 (4.41) 
and q can be obtained Be - i 
as (BY) y (4.42) 


The problem of finding (Bt) 7+ iS not a trivial 


one ana for many types of manipulators it does not exist. 
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One way of checking for the existence of the inverse is to 
compute the determinant of B; 1f its value envers esi oimiee 
is better to try to avoid it because a tremendously high con- 
trol gain will be required. 

By far the best situation is when the distributions of 
the manipulators coincide with the eigenfunctions, in which 
case B will be the identity matrix] This sreuation is mot 
very realistic; however, it is frequently possible to move the 
position of the manipulators. In this case, the best obtain- 
able positioning (n;) turns out to be the one which maximizes 
mene projection of H. (x-n,) on u, (x). 

The adequate choice of the manipulators is also a 
question of good sense, and before starting to move them and 
trying to maximize the above mentioned projection, it is high- 
ly desirable to look physically at the distribution and see if 
there 1S any possibility of obtaining it with the available 
Manipulators. As an example, it may be said that given a sys- 
tem described by Eq. 2.7 1t would be absurd to try to get the 
heat configuration shown in Fig. 4.7 with the heat source 
positioned as shown. 

It is now possible to represent the control of the 
given system in an easily implementable way, as illustrated 
iad. 4.8. 

It should be mentioned that the coefficients Ts are 
not the Fourier coefficients of the desired output distribu- 


ir 
(tafon Bale = f DS EN OCS, but these values multiplied 
0 


oF 
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aa ew we eee ee ewe = 


Figure 4.7. Example of an Undesired Control 


feed A .G. ab 
by eC. , the inverse of the closed-loop gain; 
ae) 


: represents the gain of the ate manipulator. This statement 
is also not alwaysS exact as will be seen later in the analysis 
of the computer results. Meanwhile it may be added that in 
Many cases a minimum square error deviation may not be de- 
eaeeea but another kind of criteria; in such conditions the 
t's must be chosen uSing an optimizing program, as for example 
the gradient search. When the number of filters is too high 
this search probably will require a great deal of computer 
time and it may become difficult to determine if the obtained 
optimum is local or absolute. This procedure was implemented 


successfully for three reference coefficients and one hundred 


2 Enos EOCuct Of the coectticients Tt. by the mumverse oF 


closed-loop gain will be mentioned throughout the rest of this 
thesis as the "compensated Fourier coefficients." 
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solutions of the system only took two minutes of computer 

time. 
a5 Practical) Lani ciedems 

As in the lumped parameter case the most common situa- 
tion takes place when the order of the approximated system 
differs from the number of sensors or Manipulators. In this 
case the problem of inverting non-Square matrices must be con- 
Sidered and it can be solved using the concept of the pseudo- 
inverse (¥). 

The rules for calculation of the pseudo inverse of real 
matrices are simply synthesized as follows. oie saecetiven IgE 
formation see Refs. 42 and 73. 

Ae Lt A tS. in Gra gonal fomn, AY iS also diagonal, witens ene 
non-zero elements the.reciprocal of those of A and watt 
zeros where A has zeros. 

b. Knowing that AYA 1S a square matrix obtain AY as 


av Gn A nay 


~ 


(4.43) 
The computer implementation of these concepts is shown in 

the next chapter. Depending on the numerical values of the 

elements in the matrices sometimes one procedure is better 


than the other one. ? 
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V.- DETAILED SOLUTION OF A PROBLEM BY MODAL CONTROL TECHNIQUES 


i. PROBLEM DESCRIPTION 

It is desired to heat a given rod at three different 
points such that a certain temperature distribution is 
achieved. 

Geometrically, Fig. 5.1 shows the system Con figirmotreonn 
Looking at one differential volume element it is possible to 


derive the equations describing the process. 


Heat Input (q,&t)) 


Xa arn 


Figure 5.1. The Heating of a Rod 


q(x, t ) 


SPeCtve Gaull le 

heat flow 

- cross section 

- lateral element 
of area 

av - volume element 

Op lend tomer 


OFEg 


NpaQa< 
1 








5 OG) a= 
Figure 5.2. Heat Flux in a Volume Element 


@ensider the principle of conservation of energy; from 
Fig. 5.2 the following equation can be written: 


(q,-d5)A + q3d5 = Be, (Gea) 


IL: 








However, 


qj (x4) 
and this gives 
0d) 
ee 29 — Opn (S23) 


Taking B= cpy, where c is the specific heat at constant 


volume and p the mass density it follows 





244 ayy 
= = + G39 Ax 4 =e CG qe dx (5-4) 
Or 
mee (-k ie R°dx, + 2m7Rdx, = ei Rd (oy) 
ax ee | elie ena ames 1s Sh 
it 1 1 
where k is the thermal conductivity (Btu/sec-ft-°F). Elimi- 


nating the common terms and taking a = k/pc (qa is the thermal 
rr rmsLViLty) : 


re) Yi 24. oY4 


Cm pcR TE CS 


il 
This equation can be normalized and simulated, taking into 


account the scaling factors, as 


2 
9Y tq = 2, (0<xel, €>0) Gay) 
A 2 —3t . 


For simplicity the following initial and boundary condi- 


mEvons are considered: 


Il 


y (x,0) Oy axl (By. 8), 


0 
} 0K<t <0 (Seo5 


sX(1,t) = 7 


y(0,t) 


| 
oO 
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Assume the availability of two ditferent setsworetice] 
manipulators, one giving a heat distribution coincident wren 
Ele eigenfunctions, as shown in Fig. 5.3, and the other giving 


ene heat distribution 1llustrated sine cas a: 





alemuae 1S) pS fe llacha (Dharm owiciemy Me. I 





AX 


Pique 5.4. Heat DistributaonisNo. 2 


mavoughnout the simulation runs only the two desired tem- 
perature distributions of Figs. 5.5 and 5.6 will be considered, 
the first one of which has the shape of the first eigen- 


function and being for this reason easily attainable. 


1S 





fo ee ae 


Temp 





- 4 


Figure 5.55. Temperature Distuipmrien Now 
Temp 
Z 
O 
1 x 
Figure 5.6. Temperature Distribution No. 2 


B. DERIVATION OF THE EIGENFUNCTIONS AND EIGENVALUES 
Consider again Fig. B.1l. By Laplace transforming the 


distance dependent operator with respect to t it follows that: 
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Figure B.1. Representation of the System in Operator Form 


LZ 


[--S + s]{-+} 
_ ee ae dx 
M{-] = L [-] = wis) 2 O- ee oo (5.2.10) 


To find the eigenfunctions make 


feo (Seales 
Z 
- a (s-A) u=0 
= Bre (ea 
du (xs) 
u(0,s)=0, Gee = 0 
=I) 


whose solution 1S 


jJ¥A-s -—jVA-SX 


u= Ke + Kye (55 hs) 


From the boundary condition u(0,s)=0 


Wise S) = GC Siimy «eS ey IL) 


and, from the other boundary condition, 


a TT 
dn Ss = (2n+1) 5 (Gaya lis) 
Or 
x 21° 
Be =s + (2n+1) a (DG) 
In order to normalize the eigenfunctions take pa Bem aCe oa Lor 
ay) 
f u_ (x,s) dx = 1. The following expression results: 
0 
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es [ sin? a= Xd) — a (S27) 
from which C_ = V2 
The next step 1s the determination of the adjoint operator, 
mM’, by uSing the common technique indicated in Ref. 45, pp. 269- 
feeand by Crandall [Ref. 12, p. 211]. 
By definition of adjoint operator write 
<M u, v> = <u, M’ y> (5.4.6) 


Integrating by parts the left-hand side: 


ul ge ui 1 


a we Se. du dv 
dx 0 0 
(5.19) 
LE 1 1 ib Z I 
fyemcs—-’)-vdx = [- wv] + Uy =| = | u.-2¥ + f (s-d).vdx 
0 0 Oe eb 0 


From the equations just above the definition of adjoint 


operator is verified taking 


2 
Mn = fe = re (s-) ) CS: 22.07) 
ax 
and 
slo) 2 SY 26 (5.21) 
Cite - 
x=0 


Because Mt = M and the boundary conditions for Me are the 
same as those for M, the former is said to be self-adjoint. 
This fact permits writing: 
co 1 
ee) aC) vee) cle (522) 
n=l 0 
ey accOrding to an important property of the functions of 


fmeoperator (see also Eq. Bll). 
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ee 1 
chil s)) = giggle) Gale v(e) [eld ¢ (5.23) 


which means that if £(M[-}) = M-[-] = L[-] 


mee tOllows that 
1 ii 
ayn (*) J “e) led cle (5.24) 


n 


1B] |e 


ijt=—1 8 


n=1 


Next, compute the first six eigenvalues from Eq. 5.17: 


1 Als VA. -S 
1 i 
0 S+2.4674 P25 7030 


JE s+22.2066 ee aS, 

Z Sao Oa S0 128052 8 (Say 
3 Sh L210) 5 GO) TIE) a SIS, 6 

«| So 35 Ss lee els ey 


5 S29 38RD o ec ore 


The corresponding first three eigenvalues of the operator 


L are 
- al 
NWS) So Mlerenr van 
IL 
4, (Ss) = 373272066 an 
= is 
ho(s) = sei 76850 


See oieUTATION OF THE FEEDBACK CONTROL LAW g 
Suppose it is desired to make the first three eigenvalues 


GOincide with the fourth one. Using Eq. 4.30 it follows: 


v = s + 120.9027 -s - 2.4674 = 118. 4353 
rT = sg + 120.9027 - s - 22.2066 = 98.6961 (5.27 
>, = s + 120.9027 - s - 61.6850 = 59.2177 
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D. COMPUTATION OF THE SENSORS' POSITIONS AND OF u~- 


Consider now the cases of having three and six sensors. 


The second of Eqs. 4.35 is repeated 


Oe ee 
WIGS rae eT (Cog) hh at (oon) 
EBD en) Sw AO) 
a : 
fa 3? Uy (Xe) --=- EN) ae 
where XprXg ree 1 %o are the sensors! locations. 


Two computer programs were written for the conditions 
S?N. The first one divides the interval x(0,1) in ten parts 
and tries all the possible combinations of XXX within 


this interval, such that O<X, <x, <x 3<1. If the temperature 

at x=0 was not known, the possibility of a sensor there should 
also be considered. Once, the optimal rough positions are com- 
muieea, this information is introduced in a gradient subroutine 
which searches for the rigorous optimal positions. Both pro- 
grams make use of the same function subprogram (EVMAX), which 
in turn uses two IBM library subroutines (MINV and MPRD, 
respectively for the inverse and product of matrices) and one 
N.P.G.S.'s subroutine (JACVAT) for the computation of the 
eigenvalues of a real-symmetric matrix. 


The optimal positions and corresponding U and u" matrices 


@emeeved for the cases of three and six sensors are: 


CASE I— 
x(1) to x(3) = .28570E0O0 .57144E00 205/ LEGO 


(5.258) 
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U,, to U,, = .61357E00 .13787E01 .11058E01 
U,, to U,, = .11057E01 .61353E00 -.13788E01 cao 
U,, to U,, = .13788E01 -.11057E01 .61367E00 
U,; to U,; = .17532E00 .31593E00 39391800 
U,; to U,4 = .39393E00 .17533E00 -.31591E00 (5.30) 
U,, to U,, = .31589E00 -.39391E00  .17532E00 
CASE II- 
x(1) to x(3) = .15406E00 .30828E00 .46234E00 
x(4) to x(6) = .61621E00 .76996E00 .92325E00 ae 
U,, to U,, = .338904E00 .938860E00 .132315E01 
U,, to U,, = .658373E00 .140436E01  .932899E00 
U,, to U,, = .939128E00 .116083E01 -.665084E00 (5.32) 
U,, to U,, = -116491E01 .333099E00 -.140276E01 
U,, to U,, = .132288E01 ~-.661515E00 ~-.330576E00 
Uz, to U,, = .140394E01 -.132272E01 .116496E01 
BBE to Uy, = -05225 .10148 .14464 .17915 .20310 .21531 
U5, to Uj, = .14483 .21645 .17871 .05125 -.10151 -.20289 
Us, to Ul, = .20383 .14385 -.10216 .21588 -.05109 .17894 
(S2ee) 
FE. COMPUTATION OF THE REFERENCE COEFFICIENTS 
Using the definition of fs write 
i 
tT. (t) = , i (x,t) v(x) dx (5.34) 


where i(x,t) is the desired temperature distribution. This 
gives, for the two desired temperature distributions shown in 


feos. 5.5 and 5.6, respectively 


Jyle 








TAU (1) = oe) = 0 
TAU (2) = 00.210 
TAU AGS) —— 0750 
and 
TAU (1) S35 02960 17 
TAU (2) 55 = S026 85os% 
TAU (3)) 3="-0- 2200s 
Pee DERIVATION OF B FROM THE 


~s 


MANIPULATORS 
Hei the heat ALStribution Colnci dent 


tions, B is the identity matrix. 


HEAT Dist Re uLiron 


Og ale) 7 
Sis Sa) 


(S536) 


OF THE 


with the eigenfunc- 


In the case of the heat distribution No. 2, obtain from 
Pope 4.368 (using the orthonormality properties) and Fig. 5.4 
the coefficients of B as follows: 

ORS - 
es 2 J sin 21x sinsdx = .169765 
Oro a 
b = 2 f sin 2nx sin=;dx = .363783 
2 0 2 
0.5 5 
b = #2 f{ sin 20x sin=sdx = .353678 
153 0 2 
: TT 
bo, = V2 J sin ax singdx = .600211 (5237) 
j 3 
b ny | sin 7x sin-sdx = .360126 
22 0 2 
: 5 
b = /2 f sin «x sin>sax = -.085744 
23 0 2 
1 
b,, = -f/2 J cos mx sin 5dx = .424413 
b 3 
Sy i eos crn coo en 72 
0 
. 5 
b = -/2 f cos 7x sin dx = .151576 
ce 0 2 


eleZ 








After this 0 = (ie must be computed. The result came 


~~ ~w 


out to be the following one: 


Onell) tO O:(le, 3): = la 2G - 71648 2.14944 
ere, 1) to QO(2,3): 0 0A OOS ail 2 326 (5 338) 
ie) to 0O(3,3): -89194 pe eee isi G) 2 - 88434 


eae NUMERICAL SOLUTION OF THE PARABOLEC DIFFERENTIAL EQUATION 
fiAd DESCRIBES THE SYSTEM 


The system was Simulated uSing the Crank-Nicolson method. 
The following system of equations resulted, where Jl stands 
for the instant of time J+l. The normalized distance interval 


was divided in 20 parts and the time interval was made equal 


mont /400. 
Be eo cael 17 3, 4 La “3,3 42,3 
=|. Gah 43°51 0 Moy ae Sy 
{ | { 
is es “+ ate ZOE ] 
a ade) 20) lee Digna WoT, 490,45 
| | 
=3 4ilu.. - 0 bu , 0 / | 1 
\| aay | ih eae Foe: 
wes sy) 


This system 1S solved in double precision by IBM sub- 
fomeime “DGELB." In order to take in account the fact that 
the sensors are not generally in the exact subdivisions of 
the distance domain, also IBM subroutine "ALI" was used, which 


interpolates for the sensors' positions. 


deal bys 








in Eq. 5.39 the letter uw stands for temperatune (contre — 
ily to the rest of the chapter where the temperature is 
represented by y) in order to be consistent with the computer 


program. 


H. SIMULATION RESULTS BASED IN KNOWN THEORETICAL FACTS 

Ten different simulations were performed. These simula- 
tions correspond to changes in the following parameters: 
feedback gains, Manipulators output shape and gain, reference 
Signal and number of sensors. The runs involving the use of 
strong gain saturation and a control different from the eigen- 
function control gave origin to the development of new prac- 
tical procedures. Due to the relevancy of such a fact, the 
next chapter will be entirely devoted to it. 

This thesis shows, so far as it is known, the first closed- 
loop simulation of the modal control of a distributed parameter 
system by non-optimal control techniques. In order to illus- 
trate the characteristics of this model the present section 1s 
dedicated to emphasize the results of changes in the three 
most representative parameters, namely the feedback gains, 
the direct gains and the reference coefficients. 

With exception of run No. 7 all the other ones were made 
uSing the Crank-Nicolson method, as described in Chapter I1.C. 
1. This is a brute force method that has the advantage that 
it is possible to make it quite accurate only by reducing the 
meme and distance intervals, but, in turn, it requires a 
sophisticated computer program and, consequently, a long com- 


puter time. The program written makes use of I.B.M. subroutine 
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"DGELB" for the double precision solution of a system of 
linear equations, and "ALI" which takes the temperature at 
the descritized points in the rod and interpolates these val- 
ueS in Order to get the temperatures at the sensors' positions. 
At the end of this thesis is shown the computer program used 
in run No. 10, which serves as reference program for the other 
feuns . 
Run No. 7 was done using the Bubnov-Galerkin transformation 
and this permitted writing a fairly simple program uSing a 
prediction-correction numerical solution of ordinary differen- 
tial equations [Milne, Ref. 44]. This technique can very 
easily be applied to a large variety of systems, as described 
by Foster [18]. However, when dealing with the modal control 
approach, the matrix B of the Bubnov-Galerkin method cancels 
Eae Matrix Q=B7> in the feedback loop. This makes the trans- 
formation only applicable in the case of having the manipulator 
distribution coincident with the eigenfunctions, for which B is 
the identity matrix. For this reason Foster uses the method 
only in an optimal control way, applied to the linear regu- 
lator problem, such that the mentioned transformation does 
not take place. 
1. Run No. 1: Feedback Gain Equal to Unity 

The differences from the reference program are: 

a. Dimension - U(21,200) 

b. Data - NT = 200 


er Oy OF 


e. TAU(1) 


ely PHI = Poe (3) = 8 


is: 








ee DOS 3S Lb an 


EVL(ID) = 1.57080*H*(ID-1) 
EV2(ID) = 4.71239*H*(ID-=1) 

EN 3(2D) = 70s 3S ie el ele) 
H1(ID) = 1.414214*SIN(EV1(ID)) 
H2(1D) = 1. 414204*cin (ey 2 Crp 
38H3 (ID) = 1.414214*SIN (EV3 (ID) ) 


f. First I0 1s 4, second £0 is 40 and the laste: cy uo 
Q 1s the identity matrix 
he Lhe loop *DO 41°...) was not aneluded= and sana 
allowed negative control. 
In this run it was intended to obtain the temperature 
Seem rbution No. 1 (Fig. 5.5) using manipulator control No. 1 
(Fig. 5.3); the reference, temperature coefficients were not 
divided by the respective closed-loop gains. As it can be 
Observed the output never approached the desired distribution, 
converging very slowly toward an amplitude much lower than the 
desired one; this illustrates the need for uSing the compen- 
Sated Fourier coefficients as described in Chapter IV.3. 


2. Run No. 2: Increasing the Feedback Gain of the First 
Eigenfunction 


Same statements aS in previous run except for PHI(l) = 
5.0. The convergence although faster than in the preceding 
case was still too slow. Notice that the closed-loop gain 
corresponding to the first eigenfunction decreased. This run 
illustrates how the total gain of the system decreases when 


the feedback 1S increased. 
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3. Rum No. 3: Increase of the Mantoulatomse Gans 


The differences relative to case New) are secu ne lelec ee 


Cee Are) 3.62272 


Gee Ht Gl} Sd) 


f. First I0 is 1, Second 10 is 2 and the thirasonc 
40. 


Also in statement 43 the total value was multiplied by 
20, which is equivalent to an increase of the manipulator's 
gains by the same factor. The response was very fast and at 
T=.02 secs the system was practically stabilized at the de- 
Slired temperature distribution. 

The indicated value of TAU(1) was found after having 
divided the original TAU(1l) by the respective closed-loop 
Sern. The result of this division is 3.62277, practically 
equal to the 3.62272 that gave the exact desired heat distri- 
bution. The reason for this small difference is due to the 
mteeethat the matrices P and Q are approximate trans forma- 
tions (in the specific case of eigenfunction control consider 
only the effects of P) and maybe due to round-off errors. In 
general, if N is high the empirical and theoretical values must 
be very close to each other. 


4. Run No. 4: Making the First Three Eigenvalues to 
Coincide with the Fourth One 


The statements below are the only changes in the last 


fo 
PHI(1) = 118.5353 
BIZ) ="98 6961 
Prieees) = 59.2077 
MAU(1L) = 865.502 (derived value was 85-503) 
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Figure 5.9. Run No. 3 
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Again statement 43 has the normal gain. The results were 
practically the same as in the case just considered, with the 
advantage that much less heat power was required. When trying 
to increase the gain of the manipulators by a factor of 20 the 
system became unstable. 

Some other modifications were done, as for example 
PHI (2)=PHI(3)=0 but no reasonable change was observed, which 
shows the dominant influence of the first eigenfunction in 
this specific problem; it is not so in many other cases. An- 
Other change was in the position of the sensors; for x(1)= 
ms (2)—.4 and x(3)=.5 a different corresponding matrix P 
with much greater values was obtained. AS expected, the 
means tOrmation carried on by P gave the necessary eigenfunc- 
tion information and the results are practically the same as 
when the sensors were in the optimal positions (they differ 
only in the fourth decimal place). The important difference 
between these cases resides in the fact that the optimal posi- 
tions are the ones for which the errors in the measurements 
are minimized. 


bee Rum No. 5: A different Output Distribution and 
Nonlinear Control 


In order to have temperature distribution No. 2 (Fig. 
5.6) the last run was modified as follows: 

90 DO 91 ITH=1,N1 

peer = (iH 1) 4h 

B=6.263185*xX 1 (TH) 

ern) = 1.0 - COs) 


ot CONTINUE 


Mae 
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Figure 5 .ll.- Run No. 5 








A saturated amplifier with -20<a<20 was used. It re- 
sulted that although this control is quite Saturated sam goo0 
output (at steady-state it equals the one obtained, with linear 


control) was obtained with the analytically derived t's, 


namely: 
TAU C1) Sao. 0 
DAU LZ) =. 282295 
TAU CS) ee O87 0 


Statement 40 was substituted by the segeees. 
40 DO 41 I=1,MP 
IF (ALPHA(IL) .GT.20.0) ALPHA(l)=20.0 
41 tt CALPHA(L) SLT .=2030 ) ALPHA (a) 2076 
6. Run No. 6: Larger Number of Sensors 
The number of sensors was increased from three to Six. 
This required, relative to the preceding case, the following 
changes in the computer program: 
DEMENS ION -— x(6) 72 (376), DeMe WG) 
DATA — ho—6 
xi). tO X(6) aSwan Bote > Jie ancme eas i jEic poe 
TAU (1)=116.10 
AUR 2} = 8229 3 
TAU (3)=-38. 70 
The loop 40 DO 41 _etc., was removed and again 
Statement 40 became the same as in the reference program. 
This run was done in order to illustrate the technique 


of using the pseudo-inverse matrix. 
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7s Run No. 7: The Bubnov—Galerkin Trans fommacven 
The computer program for this case is also shown in 
the es of the thesis. The loop elements are exactly as in 
paragraph five but the nonlinearity was removed. Practical- 
ly no difference is noticed in steady state, which is far 
better than expected and shows that in the present example 


mfemargner order modes of the system can be neglected. 
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Figure 5.13. Run No. 7 
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‘aVI. SIMULATION RESULTS BASED ON AGNEW Teenie ue 


It is a mathematical fact that the Fourier coefficients 
of an N terms expansion minimize the square error of the ap- 
proximation. When the matrices P and Q are high order 
matrices (theoretically infinite order) the use of the compen- 
Sated Fourier coefficients for input reference gives exactly 
the desired output. This may be the reason why Gould states 
[Ref. 26, p. 242] that typically the number of Sensors is of 
the order of fifty and the number of manipulators of the order 
of ten. When the analysis described by P is already relatively 
accurate (this happens in the example given in this thesis even 
with only a third order matrix) and the synthesis done by Q is 
also precise (1S exact for eigenfunction control) the compen- 
Sated Fourier coefficients also give origin to an output 
distribution close to the desired one. For a better under- 
Standing of what follows the above examples are nemed as ideal 
cases. However, when the accuracy of the transformations 
defined by P and Q is inadequate it is possible to improve it 
by changing the reference coefficients under the control of a 
Seacdient search until an output distribution that better 
approaches the desired one 1s achieved. In order to be able 
to use the gradient search, similarly to what was done for 
the computation of the optimal sensors' positions, a certain 


cost-function (SQUERR) was defined. 


Zo 





If it is intended to minimize the square error between 
the desired and obtained distributions, one econvenene de fini- 
Bilton Of SOQUERR iS, acconding stout q-. S alee 
M 


SQUERR = ) (yl, - y2,) 
i=l 


2 


The computer program used to implement such a technique is 


heovm in the end of the thesis. 


Temp 
a Wie 
val 
0 4 x 
pagure 6.1. How to Define the Cost Function 


For the ideal case (Q=I and P an accurate transformation) 
SQUERR was found to be .0966 and for the non-ideal one worked 
in run No. 10 (Q4I and same P as before) it was computed as 
1117. If the compensated Fourier coefficients are used 
instead of the gradient Search in any of the given non-ideal 
cases, the cost becomes much higher (1.1 in run No. 10). In 


all the examples worked out M was taken as twenty-one. 
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Another important feature of the new technique is that it 
permits, through a convenient weighting of the different terms 
in the functional cost, the desired output to be fitted in 
certain zones even more exactly than in the ideal case. One 
example of this situation is given in the run No. 9 for which 
SQUERR was defined as 


M 
SQUERR = ) 


di 
: kK, (yl, y2.) 


L 


where M is as before and 


Kk. = 1 for 2=—10toe7 wancets seems t 
k, = 10 for i=8 to 10 and 12 to 14 
k, = 100 for i=11 


1. Run No. 8 Saturated Amplifier with 0<0540 


Up to now the signal of the flux given by the manipu- 
lators was not restricted. This does not mean that the flux 
may be negative (removal of heat) but that a change of variable 
was done in the equation describing the process, such that the 
Simulated heat flow corresponds physically to a higher level 
heat flow. In some circumstances the mentioned transforma- 
tion may be undesirable and therefore it becomes necessary to 
take into account in the simulation the fact that the .control 
cannot go negative. The output was far from the desired one 
and then it was thought that maybe another set of t's would 
be better in nonlinear cases. 

To verify the suppositions just mentioned and knowing 
that the Fourier coefficients give a minimum square error 


distribution, a cost function was defined (SQUERR) which was 
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Figure 6.2. Run No. 8 








the summation of the squares of the differences between the 
desired and the obtained distributions at every one-twentieth 
Fraction of the system's dimension. The new t's did not come 
out very different but the deviation suggested that in some 
cases it 1S possible to have a better output with a set of 
fjcwoitferent from the Fourier coefficients. This is true, 
as it will be seen in the next runs. 


The optimal set of reference coefficients is 


PAU) w== sieibAr G2 
TAU(2) = 84.38 
TAU(3) = 38.70 


iiemremaining of the program is as in run No. 5, except for 
the nonlinearity which is defined as 
40 DO 41 I=1,MP 
IF (ALPHA(I) .GT.40.0) ALPHA(I)=40.0 
Al LF (ALPHA(T) 22.020) ALPHA(I)= 0.0 
meee Run No. 9: A Locally Better Output Distribution 
The nonlinearity was removed and it was decided to 
reduce the deviation between the theoretic and real outputs 
in the central zone. In order to achieve this result a 
greater weight was given to the differences between positions 
eight and fourteen, namely ten times for all except the central 
one (number eleven) which received a weight of one hundred. As 
a result the square error is not minimized but the distribution 


approaches the desired one. 
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The resultant set of T'S was: 


TAU(1) = 116.9 
TAU(2) = 18 a MA 
TAU (3) F300 


The consequences of these results may be Gulte ampeor— 
tant. Although in the present problem three eigenfunctions 
approximate already closely the required output, in other 
circumstances this may not happen. Then, the possibility of 
being able to get a closer output in some regions may indeed 
be relevant. 

Be Rum No. 10: A different Manipulator Contre! 

The computer program for this case is the one given as 
reference. As mentioned in the beginning cf this section this 
is a Situation for which the use of the compensated coeffi- 
cients is not recommended. As it can be seen in the program 
the coefficients obtained uSing the gradient search came out 


substantially different from the theoretical ones. 
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Figure 6.4, Run No. 10 
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VII. CONCLUDING REMARKS 


Aw SUGGESTIONS FOR FUTURE RESEARCH 
Many more difficult problems than the one worked out may 
appear, such as cases involving: 
1) Non-homogeneous boundary conditions; a transformation 
of variable is very useful in cases like these. 
mm DiLscontinuities. 
iii) Time varying coefficients. 
1v) Hyperbolic differential equations. [In this case, 1f 
using the Bubnov-Galerkin transformation the set of 
State variables must be augmented with those correspon- 
dent to the second derivative. Such procedure implies 
the need for knowing the initial conditions of the first 
derivatives, which is not a normal situation, and there- 
fore the use of an observer has to be implemented. 

For the beginner that gets involved in difficult situa- 
tions it may be necessary to recur to the applied mathemati- 
Cian. However, before that, there are three references that 
must be considered very carefully because of the tremendous 
potential of knowledge they contain: Courant and Hilbert 
[11], Crandall [12] and Murray-Lasso [45]. Also Foster [18], 
working in an optimal control fashion, develops a systematic 
procedure for the modeling of linear regulator problems, 


using the Bubnov-Galerkin transformation. 
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As areas for future research, problems involving non- 
linearities in the control, time variant coefficients and 
parameter identification are recommended, as well as a gener- 
alization of the Foster optimal controle method) tor some aoe 


problems other than the linear regulator. 


See CONCLUSIONS 

The most significant of the original contributions con- 
tained in the present work is, certainly, the first closed- 
seo Simulation of the modal control of a distributed 
parameter system by non-optimal control techniques. The con- 
clusions may be summarized as follows: 

1) The classical control methods can be applied to each 
independent loop. A constant in the feedback is general- 
ly good enough to obtain faster response. If it is 
desired to work close to the stability limit, design 
procedures such as Bode diagrams or root-locus plots 
also seem to be very suitable to this purpose. 

ii) The available theory states that whenever the control is 
linear (at least in steady-state), the Fourier Coeffi- 
cients of the desired output distribution, multiplied 


Gen. Ge 
by tne Laccor soon aes , Give the compensated Fourier 


X.G. 
coefficients bdew need as forcing functions (t) insure 
that the mean square error of the output is minimized. 
As shown in the computer results the above state- 
ment was proven to be not always true, specifically when 


the transformations © or P are not acecllrate, whlen 2oucme 


most general situation. Notice however that Q is an 


~~ 


INS) 





a1) 


iv) 


exact transformation (identity matrix) Sten ernescc oman 


eigenfunction contrel. Also P and @ may bevaqu2 towne = 


curate when they are low order matrices and it is neces- 


Sary to use a large number of eigenfunctions to 
approximate with precision the output or the control. 
When this happens the compensated coefficients may be 
quite far from optimal; such situation can be partially 
avoided using the technique of gradient search for the 
reference coefficients T, as originally developed in this 
thesis. 

Through an adequate definition of a cost function it is 
possible to vary the reference coefficients in such a 
way that the minimum square error is always achieved or 
even such that some regions fit the required distribu- 
tion more closely than that minimum square error 
approximation. 

The best set of manipulators is the one that gives a 
heat distribution coincident with the eigenfunctions. 

In this case G becomes the identity matrix and the sys- 
tem can be easily modeled using the Bubnov-Galerkin 
method. In many circumstances it will be difficult to 
have that type of manipulators; then, the set that ap- 
proximates more closely that distribution is the desired 


one. 
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APPENDIX A 


TRANSFORMATION OF THE PARTIAL DIFFERENTIAL EQUATION INTO AN 
ORDINARY DIFFERENTIAL EQUATION BY THE BUBNOV-GALERKIN METHOD 


The Bubnov-Galerkin method much used by Foster [Ref. 17], 
and described in great detail by Mikhlin and Smolitskiy [Ref. 
43], 1S an approximation technique for the solution of the 
equation N y-q=0, where N is a differential operator (<2 - Ly 
in this case) that must obey certain conditions. The condi- 
tions are boundedness and existence of a completely continuous 
inverse, and they will be satisfied if y(x,t) is unique and 
the coefficients of N, as well as their first derivatives are 
continuous. 

Although apparently limited, the technique applies to a 
large variety of physical situations and it seems to start 
Peeing a very important role in the control of D.P.S.‘s. It 
is a generalized and systematic way of getting the uncoupling 
effect obtained in Chapter IV. 

Start by writing an approximate solution of the partial 
differential equation in terms of the eigenfunctions: 

oo N 


18 
y(x,t) = } w(t) u(x) = ) w (t) u, (x) (At) 
n=1 n=1 


1g In the simulation run No. 7 vector w 1s represented 
praevector 2. 


? 
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Applying the condition that Ny - q is orthogonal to the 


N eigenfunctions considered here it follows that 


dw(t) 
ae = Awlt) + o(t) ae 
where 
w(t) = [w, (t) We (t) = wy (E) 1° 
e(t) = [zy(t) c(t) --- g(t)" 
(A. 3) 
$C; (t) = <q(x,t), u, (x) > 
les ye = <u; (x), L Bia Wel) 
ii. cOlwimig for C, (t) 
JL 
Eo (te) Haney (aan) u, (x) dx (A. 4) 
0 
But, 
M 
q(x,t) = ae om ft) H, (xr€,,) (A.5) 


where 0. and H were defined in Eqs. 4.56 and 4.51. Substi- 


poms A.S in A.4 it follows 


M 
2% Dad a, (t) CAuG) 


Cs (t) = 
with 
1 


es = } u, (x) H, (x, €,) dx (Aeaa 


which is the same as the previous definition of the elements 


Gm 6. 


~e 


Therefore Eq. A.2 may now be rewritten as 


dw (t) 
dt 


= Aw(t) + Ba (t) (A.8) 
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In the case of N filters and manipulators a(t) = 


qe cE a il T.- : 
(B™ ) “4H and B a(t) vie (B ) “a =U which@@eans that it is 


MOEenecesSsary to consider the mate Ea in the model. Asa 
consequence, the simulation by this method only reproduces 
the real situation when the heat distribution of the manip- 


Eacors coincides with the eigenfunctions, giving B = i 


Eee COOlVing Lor Ng 





i 92u; 
iNUe = ! u, (x) ae dix (A.9) 
and, from Chapter V, u, (x,S) = Y2 sin Y\-s x. Because M and 


Mt are self-adjoints, the normalized eigenfunction set is 
orthogonal and applying the inherent properties the matrix A 
becomes diagonal. 

A parenthesis 1s opened to say that the computation of 
the eigenfunctions done in Chapter V.B was not indispensable. 
As a matter of fact it would have been enough to choose an 
arbitrary orthogonal basis forming a complete set and satisfy- 
ing the boundary conditions, and do all the derivations with 
this basis. One set that is chosen very often and that, by 


coincidence, is the derived set of eigenfunctions is 


Zi 


B. (x) = yy Zeca 5 


(oe ~ CRS ree 





Back to the elements of A their values are: 


gael 2 
Pe, Al se TX ee es 
All = 5 5 Sin’ 5 dx = a 2.46740 
21 2 
an, = - gh fs sin? UE ax = - GP = -22.20661 
0 (Aa 
7 wed 2 
- . 2oil wnt Dix = 2 - ee 
Aj, = ae Siig) Se che = 7 = ~ 61.68503 
Ay = 0, 1 F j 
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and now all the elements in Eq. A.8 are known. The initial 
conditions are GC, (0), equal to the qe! coefficient of the ex- 
Pansion of the initial state y(x,0) andwiterstne scree ci te 


problem treated here they turn out to be zero. 


The output is given by 


y (x,t) = C(x)w(t) (A.12) 
where 


C (x) = [u, (x) uy (x) see uy Cx) J (A.13) 
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APPENDIX B 


BASIC DEFINITIONS IN FUNCTIONAL ANALYSIS 


operator 

An operator is a mapping of functions into functions. 
Figure B.1 is a block diagram representation of the operator 
fyewhiich maps the function m(x,t) into the function 7 Gee 
mine set {m} for which L is defined is the domain of L. “ine 
corresponding set iy ects range of L. An operator can 
take several configurations, the most frequent of which are 
the differential operator and its inverse, the integral 
Operator. 

When computing the inverse of an operator (m(x,t) = 
Ley (x,t)) it is necessary to know the points at which it is 
Singular and these points are called the eigenvalues hes Cr 
the direct operator (L). The eigenfunctions are then defined 


as the functions u which satisfy the equality 


Lb ° Ww =" Mei ee ee (B.1) 





Hagure B.1. Representation Of the Uncontrolled pysuen 
in Operator Form 
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The partial differential equations considered here have 
an infinite complete set of eigenfunctions and by complete, 
it 1S meant that any function with a finite number of dis- 
continuities can be represented by a linear combination of 
these eigenfunctions. This is the concept of generalized 
Fourier series. 

Peeeetnner Product 
The inner product of two continuous infinite vectors p 


and q is defined as 


b 
<prq> = f p q dx (B.2) 
a 
for finite values of the integral. If the vectors turn out 


to be discrete infinite this becomes 


_ &im 
“de a ee Pd. (B.3) 


k- 
ler S 


a 
provided that also this summation is finite. The limit is 
dropped for the case of finite vectors. 

When two sets of fectors are such that 


Oo, p7gq (B. 4) 


il 


<p ,d? 


ll 


l,p=q 
these two sets are said to be orthonormal. 
Bee j}Oint Operator 
tiie adjoint of an operator M 1s Ghe operator M’ such that 
SM u, v> = <u, M’v> (B25) 
An important property of the adjoint operator is that its 
eigenfunctions (v,) form a Set Orthogonal LO) eMercigcn Une 


tion's set Cu; ) of M. 
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If M= M* the operator is said to be Sselt-adjoimesencaean 
this case both M and M* have the same orthogonal set of 
eigenfunctions. Its properties are analogous to the proper- 
ties of the symmetric matrix and it can be proved that the 
self-adjoint case gives real eigenvalues. 

4. Spectral Representation of an Operator 
Using the orthogonality of the sets u; and Vis previously 


normalized, it follows easily that any operator with a dis- 


crete infinite spectrum may be represented as 
b co 
SoG) GYRE 9) ie) i cliche (B.6) 
a i=l 


which takes the name of spectral representation of an 


@eeractor. 


If q(x) = Lp(x), from the last equation and also from the 


foe) co 
Fourier series for p(x) = } a, u,(x) and q(x) = ) b, u, (x) 
i=} i=1 
femteuins out that 
q(x) = ee TGs) Are a. u, (x) (B.7) 
Therefore 
mea u, (x) = Pee (B. 8) 
from which, by the orthonormality properties, 
ba. = hat (ao) 
J i 


This means that the coefficients of the eigenfunction ex- 
pansion of the output can simply be Obtained LFronmenie come. 
ponding coefficients of the input by a simple product and 


therefore the operator was diagonalized. 
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A useful property of the diagonalized coeratons tose ie 


£unection of an operator whichesimp ly soeare— 


b 
ae eee 
a 


I] 1 8 


ebro rae ECE (BE 18) 


al 
For the important case of the inverse operator the spectral 
representation is 
1 2 acme il 
ST aM) OS) orem (C7 le) |i ete (B.11) 
22 mee! aL 
a i=l i 

In Ref. 45 and extensive theory of operators is given, 
which furnishes a deep knowledge about procedures for dealing 
with complex situations. This thesis will be restricted to 


separable operators for which in the case of partial differen- 


tial operators of the type 


L, _im(x,t)] 2 H.(m] + HL fm) (B.12] 


The separability is realized if 


H [HL [°]] = H, (HL 1-1) (B.13) 


Beemer restricting H, to have constant coefficients and 


G 
Laplace transforming the time, the spectral representation 
of L becomes 


L, .f'1 = } 


1 : 
oye 224 Cs) ay OO v, (¢) dz (B.14) 


A more systematic way of obtaining the diagonalization of 
an operator is using the Bubnov-Galerkin transformation, as 


described in Appendix A. 
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5. Other Defanitions 
Linear operator: an operator L is linear if the mapping 
that it implies 1s such that for arbi tramweccalancma) and 1S) 


mt Follows that Liaxs + bx.) — rN ip ae olin ca 


1 1 2 
Continuous operator: a continuous operator is character- 
ized by the fact that if a sequence of vectors tx, } converges 
to x, then the sequence of vectors {Lx} converges to Lx. 
Completely continuous operator: a linear operator is 
completely continuous if for every bounded sequence {x} in 
a linear normed space X the set {Lx} has a subsequence which 
converges in the mean to an element of X. 


Convergence in the Mean, also called strong convergence: 


the sequence {x } 1s said to converge in the mean to x when 
|| x-x,|>0, where || x|| = { : (x) dz <m (this is the case of a 
a 


12 space, which implies that x is a function defined in the 


interval a-b) or || x | =\/ ) [xf | <e (for {x} a set of complex 
== ll 


numbers and it is said that the space is an ile space). 
Bounded operator: a linear operator L is bounded if its 

domain is the entire space and if there is a scalar M such 

that || Lx| <M] x|| . The smallest of these bounds is called 


t 
norm of the operator and denoted by ; Mj . 


Subsequence: given a sequence S: {Lx, , Lx, , Lx, jee 
a subsequence of S is a fraction of it given by Lx, j Lx, 4 
i zi 
ae < < reer Vliet Sa oivane 
eens’ ee sins : where ny n. n 3 ns 
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